Diversity patterns of the Flora of Greece with applications in R
Όλα τα απαραίτητα αρχεία βρίσκονται στο github.
1 Πριν ξεκινήσουμε
Βεβαιώστε ότι το όνομα του project σας, όπως και το file path ΔΕΝ περιέχουν:
1. Ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)
Επί παραδείγματι το file path σας θα μπορούσε να είναι: C:/Environmental_Data_Analysis/Geographical_Data.
2 Φόρτωση βιβλιοθηκών
Ας φορτώσουμε όλες τις βιβλιοθήκες τις οποίες θα χρησιμοποιήσουμε σήμερα1.
## ===========================================================================##
## Load the libraries
## ===========================================================================##
library(tidyverse)
library(raster)
library(sf)
library(ape)
library(phyloregion)
library(picante)
library(rtrees)
library(ggtree)
library(exactextractr)
library(caret)
library(psych)
library(PhyloMeasures)
library(phyloregion)
library(Matrix)
library(broom)
## ===========================================================================##3 Σχέση έκτασης - αριθμού ειδών
Η σχέση αυτή θεωρείται ένας από τους ελάχιστους νόμους της οικολογίας, καθώς έχει ισχύ σε όλες τις ταξινομικές βαθμίδες και όλους τους οργανισμούς (Triantis et al., 2012). Με πολύ απλά λόγια, θεωρεί ότι όσο αυξάνεται η έκταση μιας περιοχής, τόσο περισσότερα taxa αυτή φιλοξενεί.
Όμως υπάρχουν και άλλοι παράγοντες (κλιματολογικοί, τοπογραφικοί, γεωλογικοί και εδαφικοί), οι οποίοι μπορούν να επηρεάσουν τα πρότυπα ποικιλότητας και κατανομής των ειδών στον χώρο και τον χρόνο.
4 Περιοχή μελέτης
Σήμερα θα ασχοληθούμε με τα Ryukyus Islands, τα οποία είναι περισσότερο γνωστά ως Nansei Islands (στα Ιαπωνικά: Nansei-shoto, δηλαδή “Νοτιοδυτικά Νησιά”), ανήκουν στην Ιαπωνία. Αποτελούν ένα νησιωτικό τόξο, το οποίο βρίσκεται σε γεωγραφικό πλάτος 24ο - 31.0οΝ και εκτείνεται σε μήκος 1300 km σε έναν άξονα ΝΔ-ΒΑ, ξεκινώντας από το Kyushu και καταλήγει στην Taiwan. Το νοτιότερο νησί είναι το Yonaguni, ενώ το μεγαλύτερο η Okinawa. Τα μεγαλύτερα σε μέγεθος νησιά είναι κυρίως ορεινά, ενώ τα μικρότερα είναι κοραλλιογενή. Καθώς εντοπίζεται στην ζώνη μετάβασης από το τροπικό στο θερμό εύκρατο κλίμα (δηλαδή, βρίσκεται στην υποτροπική ζώνη), οι κλιματικές συνθήκες είναι ήπιες στο αρχιπέλαγος αυτό, καθ’όλη την διάρκεια του έτους, με μέση θερμοκρασία 15οC το χειμώνα και 28οC το καλοκαίρι. Βέβαια, υπάρχουν διαβαθμίσεις: στο Β τμήμα του αρχιπελάγους, το κλίμα χαρακτηρίζεται ως υγρό υποτροπικό, ενώ στο Ν τμήμα του ως τροπικό (ταξινόμηση κατά Koppen). Εξαιτίας της γεωγραφικής τους θέσης, στα νησιά αυτά δεν υπάρχει ξηρή περίοδος και τα περισσότερα εξ αυτών δέχονται περισσότερα από 2000 m βροχής, λόγω των μουσώνων και των καλοκαιρινών τυφώνων. Ως εκ τούτου, στα Ryukyus Islands δεν εμφανίζεται η υποτροπική ξηρή ζώνη, η οποία διαχωρίζει τα εύκρατα δάση από τις υγρές τροπικές περιοχές σε άλλα μέρη του πλανήτη και η φυσική τους βλάστηση αποτελείται από πλατύφυλλα αειθαλή δαση (Good, 1845; Maekawa, 1974; Kira, 1991).
Γεωλογικά, το αρχιπέλαγος διαιρείται σε 3 τμήματα:
1. Northern Ryukyus (τα νησιά βορείως του Tokara Strait)
2. Central Ryukyus (τα νησιά μεταξύ του Tokara Strait και του Kerama Gap)
3. Southern Ryukyus (τα νησιά νοτίως του Kerama Gap)
Εικόνα 1. Γεωτεκτονική διαίρεση του αρχιπελάγους Ryukyus.
Τα Ryukyus είναι ηπειρωτικά νησιά και αποτελούν τις κορυφές της κορδιλλιέρας του αρχιπελάγους Ryukyus και έχουν δημιουργηθεί ως αποτέλεσμα της καταβύθισης της τεκτονικής πλάκας των Φιλιππίνων (Ota, 1998).
5 Εισαγωγή δεδομένων
Θα χρησιμοποιήσουμε τα ελεύθερα προσβάσιμα δεδομένα από τους Nakamura et al. (2009) για ορισμένα από τα νησιά αυτά. Τώρα ας εισαγάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε:
## ===========================================================================##
## Load our data
## ===========================================================================##
plants <- read.csv("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.csv", h = T, sep = ";",
row.names = 1)
## ===========================================================================##Ας δούμε πόσες στήλες και σειρές έχει:
## ===========================================================================##
## Check their dimensions
## ===========================================================================##
dim(plants)## [1] 26 1815
Και τι είδους δεδομένα περιέχει:
## ===========================================================================##
## Check their dimensions
## ===========================================================================##
str(plants)## 'data.frame': 26 obs. of 1815 variables:
## $ Codonacanthus_pauciflorus : int 0 0 0 1 1 0 0 1 0 0 ...
## $ Dicliptera_chinensis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Hemigraphis_reptans : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Hygrophila_salicifolia : int 0 0 0 1 0 0 0 1 0 0 ...
## $ Justicia_procumbens : int 0 0 1 1 1 1 1 1 0 1 ...
## $ Lepidagathis_formosensis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Lepidagathis_inaequalis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Staurogyne_concinnula : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Strobilanthes_glandulifer : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Strobilanthes_tashiroi : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Acer_capillipes : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Acer_insulare : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Acer_oblongum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Acer_rufinerve : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Acer_sieboldianum : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Actinidia_hypoleuca : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Actinidia_polygama : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Actinidia_rufa : int 0 0 1 1 1 1 1 1 0 0 ...
## $ Saurauria_roxburghii : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Mollugo_pentaphylla : int 0 0 1 1 1 0 0 0 0 1 ...
## $ Tetragonia_tetragonioides : int 0 1 1 1 1 1 1 1 1 1 ...
## $ Alangium_premnifolium : int 0 0 1 1 1 0 0 0 0 0 ...
## $ Alisma_canaliculatum : int 0 0 0 1 0 0 0 1 0 1 ...
## $ Caldesia_reniformis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Sagittaria_aginashi : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Sagittaria_pygmaea : int 0 0 0 1 1 0 0 0 0 0 ...
## $ Sagittaria_trifolia : int 0 0 0 1 1 1 0 0 0 0 ...
## $ Achyranthes_aspera : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Achyranthes_bidentata : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Achyranthes_longifolia : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Deeringia_polysperma : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Philoxerus_wrightii : int 0 0 0 1 1 0 0 1 0 0 ...
## $ Crinum_asiaticum : int 1 0 1 1 1 1 1 1 1 1 ...
## $ Lycoris_traubii : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ilex_crenata : int 0 0 1 1 1 0 1 1 0 0 ...
## $ Rhus_ambigua : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Rhus_javanica : int 0 1 1 1 1 1 0 0 0 0 ...
## $ Rhus_succedanea : int 0 0 1 1 1 1 0 1 0 0 ...
## $ Polyalthia_liukiuensis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Anodendron_affine : int 0 1 1 1 1 1 0 0 0 0 ...
## $ Cerbera_manghas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ecdysanthera_utilis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Neisosperma_oppositifolia : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Parsonsia_laevigata : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Trachelospermum_asiaticum : int 0 1 1 1 1 1 1 1 1 1 ...
## $ Trachelospermum_jasminoides : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Ilex_dimorphophylla : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ilex_ficoidea : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ilex_goshiensis : int 0 0 0 1 1 0 0 0 0 0 ...
## $ Ilex_hayataiana : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ilex_integra : int 1 1 1 1 1 1 1 1 1 0 ...
## $ Ilex_kusanoi : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ilex_liukiuensis : int 0 0 0 1 1 0 0 0 0 0 ...
## $ Ilex_macrocarpa : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ilex_macropoda : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Ilex_maximowicziana : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ilex_pedunculosa : int 0 0 0 1 1 0 0 0 0 0 ...
## $ Ilex_purpurea : int 0 0 0 1 0 0 0 0 0 0 ...
## $ Ilex_rotunda : int 0 0 1 1 1 1 1 1 0 0 ...
## $ Acorus_gramineus : int 0 0 1 1 1 0 0 0 0 0 ...
## $ Alocasia_cucullata : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Alocasia_gageana : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Alocasia_odora : int 0 0 1 1 1 1 1 1 0 1 ...
## $ Amorphophallus_kiusianus : int 0 0 0 1 1 0 0 1 0 0 ...
## $ Arisaema_heterocephalum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Arisaema_japonicum : int 0 0 1 1 1 0 0 1 0 0 ...
## $ Arisaema_kawashimae : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Arisaema_longipedunculatum : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Arisaema_ringens : int 0 0 1 1 1 1 1 1 1 1 ...
## $ Arisaema_robustum : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Arisaema_sazensoo : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Arisaema_thunbergii : int 0 0 1 1 1 1 1 1 0 0 ...
## $ Arisaema_yakualpinum : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Pinellia_ternata : int 0 0 1 1 1 0 0 0 0 0 ...
## $ Pinellia_tripartia : int 0 0 0 0 1 0 0 0 1 1 ...
## $ Rhaphidophora_liukiuensis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Rhaphidophora_pinnata : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Typhonium_divaricatum : int 0 0 0 1 0 0 0 0 0 0 ...
## $ Aralia_elata : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Dendropanax_trifidus : int 0 1 1 1 1 1 1 1 1 0 ...
## $ Fatsia_japonica : int 1 1 1 1 1 1 1 1 1 0 ...
## $ Hedera_rhombea : int 0 0 1 1 1 0 1 1 0 0 ...
## $ Kalopanax_pictus : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Kalopanax_septemlobus : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Schefflera_octophylla : int 0 0 1 1 1 1 1 1 0 0 ...
## $ Aristolochia_debilis : int 0 0 0 1 0 0 0 0 0 0 ...
## $ Aristolochia_kaempferi : int 0 0 0 1 1 0 0 0 0 0 ...
## $ Aristolochia_liukiuensis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Aristolochia_zollingeriana : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_caudigerum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_celsum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_dissitum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_fudsinoi : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_gelasinum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_gusuk : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_hatsushimae : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_hirsuitsepalum : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Asarum_lutchuense : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Asarum_monodoriflorum : int 0 0 0 0 0 0 0 0 0 0 ...
## [list output truncated]
6 Μετατροπή μήτρας παρουσίας/απουσίας
Όπως θα διαπιστώσατε, τα δεδομένα αυτά αποτελούν μια μήτρα παρουσίας/απουσίας ειδών (presence/absence matrix). Τέτοιου είδους δεδομένα είναι πολύ χρήσιμα για αναλύσεις σχετικά με την β-ποικιλότητα, όμως σήμερα εμείς χρειαζόμαστε μόνο τον αριθμό των ειδών σε ένα συγκεκριμένο νησί και όχι ποια είδη βρίσκονται σε κάθε νησί.
Ένας τρόπος προκειμένου να μετατρέψουμε μια μήτρα παρουσίας/απουσία ειδών στον συνολικό αριθμό ειδών τα οποία απαντώνται σε κάθε νησί, είναι ο ακόλουθος:
## ===========================================================================##
## Convert the P/A matrix to SR
## ===========================================================================##
ryukyus <- rowSums(plants) %>% enframe() %>% dplyr::rename(Island = name, SR = value)
ryukyus| Island | SR |
|---|---|
| Takeshima | 163 |
| Ioujima | 158 |
| Kuroshima | 426 |
| Tanegashima | 839 |
| Yakushima | 1025 |
| Kuchinoerabujima | 313 |
| Kuchinoshima | 336 |
| Nakanoshima | 494 |
| Gajajima | 234 |
| Tairajima | 241 |
| Suwanosejima | 267 |
| Akusekijima | 318 |
| Kodakarajima | 153 |
| Takarajima | 413 |
| Yokoatejima | 67 |
| Amamioshima | 863 |
| Kikaijima | 446 |
| Tokunoshima | 675 |
| Okinoerabujima | 568 |
| Yoronjima | 358 |
| Okinawajima | 1011 |
| Kumejima | 597 |
| Miyakojima | 552 |
| Ishigakijima | 877 |
| Iriomotejima | 887 |
| Yonagunijima | 590 |
7 Εισαγωγή γεωχωρικών δεδομένων
Καθώς έχουμε πλέον στη μορφή που θέλουμε τα βιοτικά μας δεδομένα, επόμενο μας βήμα είναι να συγκεντρώσουμε αβιοτικά δεδομένα για τα νησιά αυτά.
7.1 Διανυσματικά δεδομένα
Για τον σκοπό αυτό, θα χρειαστεί πρώτα να εισαγάγουμε τα διανυσματικά αρχεία (shapefiles) για τα νησιά που μας ενδιαφέρουν.
Θα χρειαστεί να αλλάξετε το file path σε
Shapefiles.
## ===========================================================================##
## Create a vector containing just the file paths
## ===========================================================================##
shps <- fs::dir_ls("C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/",
glob = "*.shp")
## ===========================================================================##
## ===========================================================================##
## Create a vector containing just the island names
## ===========================================================================##
island_nms <- shps %>% str_replace_all("C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/",
"") %>% str_replace_all("\\.shp", "")
## ===========================================================================##
## ===========================================================================##
## Bulk load the shapefiles
## ===========================================================================##
ryukyus_shps <- shps %>% purrr::map(read_sf) %>% lapply(., dplyr::select, NAME_1,
geometry) %>% do.call(rbind, .) %>% mutate(Island = island_nms) %>% dplyr::select(Island,
geometry)
## ===========================================================================##7.2 Ψηφιδοπλέγματα
Ας είσαγουμε τώρα και τα ψηφιδοπλέγματα τα οποία περιέχουν τις αβιοτικές μας μεταβλητές.
Αναλόγως των δυνατοτήτων του Η/Υ σας, υπάρχει περίπτωση οι ακόλουθες εντολές να διαρκέσουν μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτών με την εντολή
preds <- readRDS('Data/RDS/Ryukyus abiotic data.rds')&altitude <- readRDS('Data/RDS/Ryukyus altitudinal data.rds')
## ===========================================================================##
## Load, crop and mask the rasters
## ===========================================================================##
preds <- stack(stack(list.files(path = "I:/SDM DATA/WORLDCLIM/Current/bio_2-5m_bil/",
pattern = ".bil", recursive = T, full.names = T)) %>% crop(., ryukyus_shps, snap = "in") %>%
mask(., ryukyus_shps), stack(list.files(path = "I:/SDM DATA/Envirem_Predictors/Eurasia_current_2.5arcmin_geotiff/",
pattern = ".tif", recursive = T, full.names = T)) %>% crop(., ryukyus_shps, snap = "in") %>%
mask(., ryukyus_shps))
## ===========================================================================##Ας μεταφορτώσουμε και υψομετρικά δεδομένα για τα νησιά αυτά:
## ===========================================================================##
## Load, crop and mask the altitudinal raster
## ===========================================================================##
altitude <- raster("I:/SDM DATA/Altitude/elevation_1KMmd_GMTEDmd.tif") %>% crop(.,
ryukyus_shps, snap = "in") %>% mask(., ryukyus_shps)
## ===========================================================================##8 Εξαγωγή γεωχωρικών δεδομένων
Τώρα μπορούμε να εξάγουμε τα αβιοτικά αυτά δεδομένα για κάθε νησί που εμπεριέχεται στην ανάλυση μας:
## ===========================================================================##
## Extract abiotic info
## ===========================================================================##
preds_ryukyus <- exact_extract(preds, ryukyus_shps, "median") %>% mutate(Island = ryukyus_shps$Island) %>%
full_join(., exact_extract(altitude, ryukyus_shps, "median") %>% enframe() %>%
mutate(Island = ryukyus_shps$Island) %>% rename(Elevation = value) %>% dplyr::select(Island,
Elevation)) %>% dplyr::select(Island, everything()) %>% as_tibble()##
|
| | 0%
|
|=================================== | 50%
|
|======================================================================| 100%
##
|
| | 0%
|
|=================================== | 50%
|
|======================================================================| 100%
## ============================================================================##
## Give them proper names
## ============================================================================##
col_names <- preds_ryukyus %>% colnames() %>% enframe() %>% mutate(value = str_replace_all(value,
"median\\.", "")) %>% mutate(value = str_replace_all(value, "current_2\\.5arcmin_",
"")) %>% pull(value)
## ===========================================================================##
## ============================================================================##
## Overwrite names and remove Row 1
## ============================================================================##
vars <- preds_ryukyus %>% set_names(col_names) %>% dplyr::select(Island, everything()) %>%
as_tibble()9 Αντιμετώπιση προβλημάτων
Εάν παρατηρήσετε τα δεδομένα για τα νησιά που μας ενδιαφέρουν, θα δείτε2, ότι για ορισμένα νησιά δεν έχουμε καθόλου (NA) δεδομένα. Αυτό μπορεί να οφείλεται είτε στο γεγονός ότι τα νησιά αυτά μπορεί να είναι πολύ μικρά, οπότε να μην υπάρχουν διαθέσιμα δεδομένα στο μέγεθος των ψηφιδοπλεγμάτων τα οποία χρησιμοποιούμε (resolution - ακρίβεια), είτε εξαιτίας του γεγονότος ότι χρησιμοποιούμε σχετικά αδρά ψηφιδοπλέγματα να μην έχουμε δεδομένα για αυτά τα νησιά, ανεξαρτήτως του μεγέθους τους3.
Υπάρχουν δύο αποδεκτές λύσεις:
- Μετασχηματισμός των δεδομένων μας
- Αφαίρεση των προβληματικών δεδομένων
Και οι δύο προσεγγίσεις έχουν τα πλεονεκτήματα και τα μειονεκτήματα τους. Η πρώτη προσέγγιση είναι προτιμότερη όταν έχουμε πλήθος διαθέσιμων δεδομένων (>5000 δείγματα) και το ποσοστό των προβληματικών δεδομένων είναι μικρό (συνήθως < 1-2%). Η δεύτερη προσέγγιση θεωρείται γενικά ασφαλέστερη, διότι δεν δημιουργούμε εκ του μηδενός δεδομένα. Σήμερα θα ακολουθήσουμε την δεύτερη προσέγγιση, όμως θα δούμε πως μπορούμε να μετασχηματίσουμε τα δεδομένα μας.
9.1 Μετασχηματισμός δεδομένων
Ας δούμε πρώτα πως μπορούμε να μετασχηματίσουμε τα δεδομένα μας:
9.2 Αφαίρεση δεδομένων
Ας δούμε τώρα πως μπορούμε να αφαιρέσουμε τα όποια προβληματικά δεδομένα μας:
## ============================================================================##
## First drop any rows with NA values
## ============================================================================##
drop_vars <- vars %>% drop_na()
## ============================================================================##
## ============================================================================##
## Then check which islands have NA values in any of the variables
## ============================================================================##
setdiff(imputed_vars$Island, drop_vars$Island)## [1] "Akusekijima" "Gajajima" "Ioujima" "Kodakarajima" "Kuchinoshima"
## [6] "Tairajima" "Takarajima" "Takeshima" "Yokoatejima" "Yonagunijima"
## ============================================================================##
## Create a symbol to filter out all the Islands having NA values
## ============================================================================##
"%!in%" <- function(x, y) !(x %in% y)
## ============================================================================##
## ============================================================================##
## Now create a temporary presence/absence matrix removing any problematic islands
## ============================================================================##
plants_new <- plants %>% mutate(Island = rownames(.)) %>% as_tibble() %>% dplyr::select(Island,
everything()) %>% pivot_longer(!Island, names_to = "Taxon", values_to = "PA") %>%
filter(Island %!in% setdiff(imputed_vars$Island, drop_vars$Island)) %>% pivot_wider(names_from = Taxon,
values_from = PA)
## ============================================================================##
## ============================================================================##
## See which taxa are left after removing the problematic islands
## ============================================================================##
keepers <- colSums(plants_new[, -1]) %>% enframe() %>% filter(value > 0)
## ============================================================================##
## ============================================================================##
## Finally create the new PA matrix
## ============================================================================##
plants_new_pa <- plants_new %>% pivot_longer(!Island, names_to = "Taxon", values_to = "PA") %>%
filter(Taxon %in% keepers$name) %>% pivot_wider(names_from = Taxon, values_from = PA)
## ============================================================================##10 Κατασκευή φυλογενετικού δένδρου
Προκειμένου να μπορέσουμε να υπολογίσουμε την φυλογενετική α-ποικιλότητα για το αρχιπέλαγος της Ιαπωνίας, θα χρειαστεί να κατασκευάσουμε πρώτα το φυλογενετικό δένδρο των ειδών που απαρτίζουν τις νησιωτικές φυτοκοινότητες. Για τον σκοπό αυτό, θα χρησιμοποιήσουμε εντολές από τις βιβλιοθήκες tidyr, dplyr & rtrees. Η τελευταία βιβλιοθήκη περιλαμβάνει τα φυλογενετικά δένδρα για αρκετούς φυτικούς οργανισμούς, καθώς και για τους ιχθύες, τα πτηνά και τα θηλαστικά. Θα χρησιμοποιήσουμε μόνο τα taxa τα οποία απαντώνται στα 16 νησιά του αρχιπελάγους της Ιαπωνίας που δεν είχαν κάποιο πρόβλημα με τις αβιοτικές μεταβλητές. Τα taxa αυτά βρίσκονται στον πίνακα που ονομάσαμε keepers. Θα χρειαστεί να προσθέσουμε στον πίνακα αυτόν το όνομα του γένους για κάθε φυτικό taxon το οποίο απαντάται στα 16 αυτά νησιά των Ryukyus και στη συνέχεια, θα χρησιμοποιήσουμε αυτόν τον πίνακα για να κατασκευάσουμε το φυλογενετικό μας δένδρο. Ακόμα και εάν το φυλογενετικό μας δένδρο περιέχει πολυτομίες, είναι ασφαλές να το χρησιμοποιήσουμε για μακροοικολογικές αναλύσεις (Li et al., 2019).
Αναλόγως των δυνατοτήτων του Η/Υ σας, υπάρχει περίπτωση η ακόλουθη εντολή να διαρκέσει μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτής με την εντολή:
jpn_tree <- readRDS('Data/RDS/16 Ryukyus tree.rds')
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
library(rtrees)
jpn_tree <- keepers %>%
separate(name, into = c('genus', 'binomial'), sep = '_', remove = F) %>%
rename(species = name) %>%
dplyr::select(species, genus) %>%
get_tree(sp_list = .,
tree = tree_plant_otl, # either
taxon = "plant", # or
scenario = "S1",
show_grafted = FALSE)
##============================================================================##Στη συνέχεια, μπορούμε να αναπαραστήσουμε γραφικά το φυλογενετικό δένδρο, δίνοντας διαφορετικό χρώμα σε κάθε γένος:
## ============================================================================##
## Plot the tree
## ============================================================================##
groupInfo <- split(jpn_tree$tip.label, gsub("_\\w+", "", jpn_tree$tip.label))
groups_jpn <- groupOTU(jpn_tree, groupInfo)
ggtree(groups_jpn, aes(color = group), layout = "circular") + geom_tiplab2(size = 0.5,
align = T)10.1 Υπολογισμός φυλογενετικής α-ποικιλότητας
Στο σημείο αυτό, θα χρησιμοποιήσουμε μια εντολή από την βιβλιοθήκη picante, προκειμένου να κρατήσουμε μόνο εκείνα τα είδη τα οποία είναι κοινά μεταξύ του φυλογενετικού δένδρου και της μήτρας παρουσίας/απουσίας:
## ============================================================================##
## Combine the trees with the P/A data
## ============================================================================##
# check for mismatches/missing species
combined_japan <- picante::match.phylo.comm(jpn_tree, plants_new_pa[, 2:1800])## [1] "Dropping taxa from the community because they are not present in the phylogeny:"
## [1] "Saurauria_roxburghii" "Carmona_retusa"
## [3] "Saionia_shinzatoi" "Crataeva_falcata"
## [5] "Pseudostellaris_heterantha" "Dispanthus_uniflorus"
## [7] "Siegesbeckia_glabrescens" "Siegesbeckia_orientalis"
## [9] "Cardamina_flexuosa" "Pogonantherum_crinitum"
## [11] "Spodipogon_cotulifer" "Eusteralis_stellata"
## [13] "Perillura_reptans" "Akebi_quinata"
## [15] "Dumbaria_villosa" "Alectrorurus_yedoensis"
## [17] "Caphalanthera_falcata" "Heterozeuxine_nervosa"
## [19] "Heterozeuxine_odorata" "Luisaerides_liukiuensis"
## [21] "Vexillabium_yakushimense" "Hydyotis_verticillata"
## [23] "phtheirospermum_japonicum" "Cnidum_japonicum"
## [1] "Dropping tips from the tree because they are not present in the community data:"
## [1] "Phtheirospermum_japonicum"
## ============================================================================##
## The resulting object is a list with $phy and $comm elements replace our
## original data with the sorted/matched data
## ============================================================================##
phy_japan <- combined_japan$phy
comm_japan <- combined_japan$comm
## ============================================================================##Τώρα μπορούμε να υπολογίσουμε την φυλογενετική α-ποικιλότητα, μέσω της βιβλιοθήκης PhyloMeasures:
## ============================================================================##
## Compute phylogenetic alpha diversity
## ============================================================================##
## For the estimation of SES values
japan_pd <- pd.query(phy_japan, comm_japan, TRUE)
## For the p-values
japan_pd_p <- pd.pvalues(phy_japan, comm_japan, reps = 1000)
## ============================================================================##11 Τελικό dataset
Τώρα πρέπει να κρατήσουμε μόνο εκείνα τα νησιά που δεν έχουν προβληματικά δεδομένα, να τα ενώσουμε με τα αβιοτικά και τα φυλογενετικά δεδομένα, καθώς και να υπολογίσουμε την έκταση του εκάστοτε νησιού:
##============================================================================##
## Create the final dataset
##============================================================================##
ryukyus_new <- ryukyus %>%
filter(Island %in% plants_new_pa$Island) %>%
full_join(., drop_vars) %>%
full_join(., ryukyus_shps) %>%
filter(Island %!in% setdiff(imputed_vars$Island,
drop_vars$Island)) %>%
st_as_sf() %>%
## Lambert Azimuthal Equal Area projection
st_transform('+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs') %>%
## Create the Area and PD variables
mutate(Area = as.numeric(units::set_units(st_area(.), km^2)),
PD = japan_pd,
PD_P = japan_pd_p) %>%
## Transform back to WGS84
st_transform(4326) %>%
dplyr::select(Island, SR, PD, PD_P, Area, Elevation, everything()) %>%
## It has zero-variance
dplyr::select(-aridityIndexThornthwaite) %>%
## Drop the geometry
st_drop_geometry() %>%
## Join these data with the Maximum and Minimum distance between each island
## pair
inner_join(.,
## Calculate the Maximum distance between each island pair
apply((st_distance(ryukyus_shps %>%
st_centroid(),
ryukyus_shps %>%
st_centroid())/1e+3) %>%
as_tibble() %>%
set_names(ryukyus_shps$Island) %>%
units::drop_units() %>%
na_if(., 0), 1, max, na.rm = T) %>%
enframe() %>%
rename(MaxD = value) %>%
mutate(Island = ryukyus_shps$Island) %>%
dplyr::select(Island, MaxD) %>%
full_join(.,
## Calculate the Minimum distance between each island
## pair
apply((st_distance(ryukyus_shps %>%
st_centroid(),
ryukyus_shps %>%
st_centroid())/1e+3) %>%
as_tibble() %>%
set_names(ryukyus_shps$Island) %>%
units::drop_units() %>%
na_if(., 0), 1, min, na.rm = T) %>%
enframe() %>%
rename(MinD = value) %>%
mutate(Island = ryukyus_shps$Island) %>%
dplyr::select(Island, MinD)))
##============================================================================##
ryukyus_new| Island | SR | PD | PD_P | Area | Elevation | bio1 | bio10 | bio11 | bio12 | bio13 | bio14 | bio15 | bio16 | bio17 | bio18 | bio19 | bio2 | bio3 | bio4 | bio5 | bio6 | bio7 | bio8 | bio9 | annualPET | climaticMoistureIndex | continentality | embergerQ | growingDegDays0 | growingDegDays5 | maxTempColdest | minTempWarmest | monthCountByTemp10 | PETColdestQuarter | PETDriestQuarter | PETseasonality | PETWarmestQuarter | PETWettestQuarter | thermicityIndex | MaxD | MinD |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Kuroshima | 426 | -1.8402810 | 0.9660340 | 16.30184 | 213.81479 | 176.0000 | 253.0000 | 102.00000 | 2636.000 | 516.0000 | 117.0000 | 48.00000 | 1079.0000 | 403.0000 | 717.0000 | 408.0000 | 68.00000 | 29.00000 | 5816.000 | 293.0000 | 59.00000 | 234.0000 | 221.0000 | 157.00000 | 1000.6500 | 0.6200000 | 16.90000 | 387.3800 | 76194.00 | 76194.00 | 125.0000 | 229.0000 | 10.000000 | 44.86000 | 58.52000 | 3070.150 | 113.6200 | 115.83000 | 360.0000 | 983.8313 | 34.73631 |
| Tanegashima | 839 | -2.0300010 | 0.9830170 | 449.44967 | 89.10355 | 189.1470 | 264.5754 | 115.14228 | 2759.873 | 492.1473 | 108.4328 | 44.06544 | 1075.9171 | 374.1039 | 805.9562 | 374.1039 | 61.59619 | 26.50000 | 5791.570 | 302.0701 | 73.81043 | 228.0324 | 234.1318 | 115.14228 | 990.7673 | 0.6451378 | 16.82070 | 417.7372 | 81746.30 | 81746.30 | 135.3028 | 242.8109 | 11.500000 | 44.90803 | 44.90803 | 3037.074 | 113.2602 | 114.93804 | 398.3696 | 1040.1992 | 51.08020 |
| Yakushima | 1025 | -0.2158943 | 0.5644356 | 508.57513 | 606.53490 | 158.5849 | 231.3554 | 86.39672 | 3318.087 | 572.3122 | 137.1453 | 41.59947 | 1270.1959 | 486.8870 | 964.3652 | 486.8870 | 61.16105 | 27.14583 | 5602.943 | 267.2355 | 46.68151 | 220.8119 | 199.3229 | 86.39672 | 909.2751 | 0.7249949 | 16.21991 | 520.0648 | 68752.66 | 68752.66 | 106.3139 | 210.1564 | 9.506935 | 40.78871 | 40.78871 | 2838.298 | 103.5724 | 105.73511 | 311.2719 | 989.5511 | 32.70489 |
| Kuchinoerabujima | 313 | -1.6912781 | 0.9610390 | 37.17430 | 158.07548 | 184.0000 | 256.0000 | 114.00000 | 2812.000 | 543.0000 | 123.0000 | 45.00000 | 1105.0000 | 459.0000 | 724.0000 | 459.0000 | 64.00000 | 29.00000 | 5478.000 | 293.0000 | 74.00000 | 219.0000 | 196.0000 | 114.00000 | 1000.3900 | 0.6400000 | 15.80000 | 440.4100 | 79416.00 | 79416.00 | 135.0000 | 232.0000 | 12.000000 | 45.54000 | 45.54000 | 3059.810 | 113.1800 | 107.40000 | 393.0000 | 975.0262 | 32.70489 |
| Nakanoshima | 494 | -1.0391167 | 0.8381618 | 35.72809 | 244.47767 | 185.0000 | 254.0000 | 118.00000 | 3042.000 | 564.0000 | 146.0000 | 43.00000 | 1199.0000 | 498.0000 | 748.0000 | 498.0000 | 61.00000 | 28.00000 | 5234.000 | 291.0000 | 80.00000 | 211.0000 | 194.0000 | 118.00000 | 978.1600 | 0.6800000 | 15.30000 | 494.1600 | 80010.00 | 80010.00 | 138.0000 | 233.0000 | 12.000000 | 45.69000 | 45.69000 | 2878.560 | 110.3400 | 103.21000 | 403.0000 | 906.5470 | 14.56124 |
| Suwanosejima | 267 | -0.9998780 | 0.8271728 | 29.85171 | 191.07586 | 191.0000 | 260.0000 | 125.00000 | 3033.000 | 547.5000 | 148.5000 | 42.00000 | 1182.0000 | 499.5000 | 743.5000 | 499.5000 | 60.00000 | 28.00000 | 5202.000 | 295.5000 | 86.50000 | 209.0000 | 199.0000 | 125.00000 | 988.7250 | 0.6750000 | 15.10000 | 496.4800 | 82503.00 | 82503.00 | 144.5000 | 237.5000 | 12.000000 | 46.77500 | 46.77500 | 2862.030 | 111.2700 | 103.63500 | 422.0000 | 880.1082 | 18.14204 |
| Amamioshima | 863 | -0.8816950 | 0.8081918 | 724.01156 | 168.94980 | 207.2345 | 269.9245 | 144.86514 | 2927.709 | 412.3831 | 145.5170 | 30.46525 | 964.3180 | 490.3122 | 837.2828 | 490.3122 | 57.82408 | 29.00000 | 4885.349 | 305.1555 | 109.35008 | 196.0141 | 219.1768 | 144.86514 | 1019.7703 | 0.6466562 | 13.94785 | 505.2497 | 89714.96 | 89714.96 | 165.8180 | 248.5766 | 12.000000 | 50.93358 | 50.93358 | 2752.129 | 112.8781 | 104.93128 | 482.6365 | 771.3206 | 53.54787 |
| Kikaijima | 446 | -2.4427093 | 0.9930070 | 59.58188 | 34.49266 | 213.5000 | 276.5000 | 153.01189 | 2868.929 | 416.6492 | 143.9881 | 31.50000 | 959.9048 | 479.0000 | 802.7758 | 479.0000 | 58.50000 | 29.00000 | 4851.000 | 312.0000 | 115.50000 | 196.5000 | 226.0119 | 153.01189 | 1050.6402 | 0.6350000 | 13.84612 | 494.0876 | 93132.43 | 93132.43 | 174.0119 | 253.5000 | 12.000000 | 52.36262 | 52.36262 | 2846.930 | 116.6501 | 108.38107 | 506.0119 | 817.8176 | 53.54787 |
| Tokunoshima | 675 | -1.0657777 | 0.8661339 | 253.48894 | 123.09407 | 213.5393 | 273.4286 | 153.53926 | 2409.559 | 319.2652 | 123.8041 | 27.50000 | 811.3153 | 411.0024 | 660.5954 | 411.0024 | 53.50000 | 28.00000 | 4683.245 | 305.9859 | 119.45537 | 186.5074 | 224.4286 | 153.53926 | 999.6211 | 0.5852746 | 13.27762 | 438.1796 | 92427.59 | 92427.59 | 173.5393 | 253.0934 | 12.000000 | 51.34734 | 51.34734 | 2598.477 | 110.1598 | 101.82562 | 506.3661 | 700.2143 | 55.69807 |
| Okinoerabujima | 568 | -1.8609689 | 0.9640360 | 95.13222 | 72.94901 | 217.8094 | 276.4325 | 159.32750 | 2148.366 | 273.9207 | 111.8200 | 26.50000 | 729.3018 | 368.1002 | 580.7622 | 370.2275 | 51.50000 | 28.00000 | 4602.111 | 307.4325 | 125.80945 | 181.5000 | 228.7475 | 159.32750 | 994.5787 | 0.5364186 | 12.92500 | 401.9761 | 94279.97 | 94279.97 | 179.1342 | 255.8094 | 12.000000 | 52.08294 | 52.08294 | 2518.234 | 109.2814 | 100.45288 | 522.4284 | 648.5441 | 40.89212 |
| Yoronjima | 358 | -2.2220711 | 0.9880120 | 20.92543 | 26.13938 | 222.0000 | 279.0000 | 163.00000 | 2139.000 | 274.0000 | 111.0000 | 28.00000 | 713.0000 | 362.0000 | 603.0000 | 362.0000 | 53.00000 | 28.00000 | 4556.000 | 311.0000 | 128.00000 | 183.0000 | 275.0000 | 163.00000 | 1026.5900 | 0.5200000 | 12.85000 | 396.0200 | 95922.00 | 95922.00 | 184.0000 | 258.0000 | 12.000000 | 54.10000 | 54.10000 | 2567.800 | 112.4200 | 115.28000 | 534.0000 | 616.5526 | 40.89212 |
| Okinawajima | 1011 | -1.5138595 | 0.9350649 | 1236.14528 | 74.98082 | 220.2019 | 276.5807 | 161.14858 | 2216.502 | 290.1055 | 115.5279 | 31.38394 | 777.4536 | 367.2347 | 665.8982 | 367.2347 | 54.85507 | 29.50000 | 4513.869 | 309.0611 | 124.09077 | 183.2771 | 273.3103 | 161.14858 | 1038.3087 | 0.5243882 | 12.79125 | 404.2533 | 95252.11 | 95252.11 | 182.3763 | 255.1475 | 12.000000 | 55.32691 | 55.32691 | 2511.235 | 111.9716 | 114.67432 | 526.7413 | 548.4580 | 76.76997 |
| Kumejima | 597 | -1.4178822 | 0.9250749 | 60.71743 | 54.01407 | 222.0000 | 277.0000 | 163.00000 | 2303.000 | 315.0000 | 132.5000 | 29.00000 | 809.5000 | 423.5000 | 576.5000 | 423.5000 | 52.00000 | 28.00000 | 4482.000 | 309.0000 | 128.50000 | 180.5000 | 260.0000 | 163.00000 | 1022.8400 | 0.5550000 | 12.77500 | 432.4150 | 95931.00 | 95931.00 | 182.0000 | 257.0000 | 12.000000 | 54.14000 | 54.14000 | 2539.940 | 111.3300 | 111.44500 | 532.5000 | 623.7271 | 119.42800 |
| Miyakojima | 552 | -2.3345104 | 0.9870130 | 165.13666 | 42.67389 | 235.5000 | 281.2234 | 184.50000 | 2144.536 | 242.3613 | 136.8159 | 20.50000 | 646.8572 | 417.0039 | 621.6921 | 429.3743 | 49.00000 | 30.00000 | 3854.703 | 314.2234 | 153.89129 | 160.5653 | 270.2234 | 185.89129 | 1041.4816 | 0.5150000 | 10.99898 | 451.9025 | 101817.17 | 101817.17 | 201.1930 | 261.2234 | 12.000000 | 56.05820 | 66.00715 | 2539.382 | 111.8939 | 97.58236 | 590.6965 | 850.5334 | 120.19214 |
| Ishigakijima | 877 | -0.6252511 | 0.7242757 | 230.48583 | 40.31216 | 236.1942 | 281.6206 | 185.25599 | 2132.099 | 236.0820 | 122.4948 | 21.94050 | 655.8379 | 388.7829 | 619.1044 | 400.0119 | 50.50000 | 31.00000 | 3833.381 | 313.2880 | 153.19419 | 160.0733 | 269.2880 | 187.28797 | 1055.7542 | 0.5071154 | 10.98712 | 447.8451 | 102097.75 | 102097.75 | 203.7891 | 263.2880 | 12.000000 | 58.52420 | 67.71221 | 2336.965 | 113.6409 | 99.85277 | 593.6196 | 954.2250 | 41.09759 |
| Iriomotejima | 887 | 0.1921938 | 0.3906094 | 292.10221 | 159.56565 | 228.9828 | 273.2398 | 178.98280 | 2370.300 | 258.4808 | 151.5683 | 18.23086 | 729.2650 | 471.4232 | 636.3513 | 479.9272 | 51.50000 | 32.12202 | 3736.810 | 305.5627 | 147.80150 | 158.2256 | 261.2398 | 183.77527 | 1050.9475 | 0.5553407 | 10.68895 | 506.5259 | 99032.11 | 99032.11 | 197.3049 | 252.8150 | 12.000000 | 57.64614 | 67.09907 | 2405.794 | 114.4487 | 99.20065 | 574.2911 | 989.3600 | 41.09759 |
12 Προκαταρτικές αναλύσεις
12.1 Έλεγχος κανονικότητας
12.1.1 Α
Δημιουργούμε ένα αρχείο το οποίο ελέγχει εάν όλες οι μεταβλητές μας (εκτός της πρώτης στήλης, η οποία είναι ποσοτική μεταβλητή - το όνομα των νησιών), έχουν κανονική κατανομή:
## ============================================================================##
## Normality check
## ============================================================================##
normality <- lapply(ryukyus_new[, 2:42], shapiro.test)
## ============================================================================##12.1.2 Β
Ας δούμε τα αποτελέσματα του προηγούμενου βήματος
## $SR
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.93472, p-value = 0.2893
##
##
## $PD
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.96192, p-value = 0.6968
##
##
## $PD_P
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.77558, p-value = 0.001314
##
##
## $Area
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.75456, p-value = 0.0007212
##
##
## $Elevation
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.72274, p-value = 0.0003039
##
##
## $bio1
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.93352, p-value = 0.2771
##
##
## $bio10
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.85281, p-value = 0.01499
##
##
## $bio11
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.93575, p-value = 0.3001
##
##
## $bio12
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.90384, p-value = 0.09261
##
##
## $bio13
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.86312, p-value = 0.02138
##
##
## $bio14
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92022, p-value = 0.1701
##
##
## $bio15
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92241, p-value = 0.1845
##
##
## $bio16
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.91524, p-value = 0.1414
##
##
## $bio17
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.88965, p-value = 0.055
##
##
## $bio18
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92892, p-value = 0.2344
##
##
## $bio19
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.89273, p-value = 0.06153
##
##
## $bio2
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.94329, p-value = 0.3913
##
##
## $bio3
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.91761, p-value = 0.1544
##
##
## $bio4
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.94292, p-value = 0.3863
##
##
## $bio5
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.82406, p-value = 0.005789
##
##
## $bio6
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.93515, p-value = 0.2938
##
##
## $bio7
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.94416, p-value = 0.4031
##
##
## $bio8
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.88975, p-value = 0.0552
##
##
## $bio9
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92902, p-value = 0.2352
##
##
## $annualPET
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.88922, p-value = 0.05414
##
##
## $climaticMoistureIndex
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92624, p-value = 0.2124
##
##
## $continentality
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.93748, p-value = 0.3193
##
##
## $embergerQ
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.91509, p-value = 0.1406
##
##
## $growingDegDays0
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92894, p-value = 0.2345
##
##
## $growingDegDays5
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92894, p-value = 0.2345
##
##
## $maxTempColdest
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92989, p-value = 0.2428
##
##
## $minTempWarmest
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.88613, p-value = 0.04841
##
##
## $monthCountByTemp10
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.46956, p-value = 1.199e-06
##
##
## $PETColdestQuarter
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9442, p-value = 0.4037
##
##
## $PETDriestQuarter
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92498, p-value = 0.2028
##
##
## $PETseasonality
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92548, p-value = 0.2065
##
##
## $PETWarmestQuarter
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.87995, p-value = 0.03873
##
##
## $PETWettestQuarter
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.91681, p-value = 0.1498
##
##
## $thermicityIndex
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.93209, p-value = 0.2631
##
##
## $MaxD
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.91265, p-value = 0.1284
##
##
## $MinD
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8222, p-value = 0.005454
Οι μεταβλητές με πιθανότητα μικρότερη από 0.05 ΔΕΝ εμφανίζουν κανονική κατανομή
12.2 Έλεγχος συσχέτισης
Ας ελέγξουμε τη συσχέτιση μεταξύ των μεταβλητών μας
Εάν κάποιες έχουν συσχέτιση \(\geq 0.7\), πρέπει να τις αφαιρέσουμε από τη μεταγενέστερη ανάλυση (Γιατί;)
## ============================================================================##
## Correlation check
## ============================================================================##
corr.test(ryukyus_new[, 5:42], method = "spearman")## Call:corr.test(x = ryukyus_new[, 5:42], method = "spearman")
## Correlation matrix
## Area Elevation bio1 bio10 bio11 bio12 bio13 bio14 bio15
## Area 1.00 0.06 0.17 0.11 0.14 -0.01 -0.22 0.01 -0.32
## bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6
## Area -0.15 -0.07 0.25 -0.05 -0.22 0.11 -0.22 0.10 0.10
## bio7 bio8 bio9 annualPET climaticMoistureIndex
## Area -0.20 0.24 0.00 0.11 -0.01
## continentality embergerQ growingDegDays0 growingDegDays5
## Area -0.22 0.32 0.16 0.16
## maxTempColdest minTempWarmest monthCountByTemp10
## Area 0.15 0.09 -0.09
## PETColdestQuarter PETDriestQuarter PETseasonality
## Area 0.18 -0.05 -0.44
## PETWarmestQuarter PETWettestQuarter thermicityIndex MaxD
## Area -0.06 -0.23 0.14 0.04
## MinD
## Area 0.49
## [ reached getOption("max.print") -- omitted 37 rows ]
## Sample Size
## [1] 16
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## Area Elevation bio1 bio10 bio11 bio12 bio13 bio14 bio15
## Area 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
## bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6 bio7
## Area 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
## bio8 bio9 annualPET climaticMoistureIndex continentality
## Area 1.00 1.00 1.00 1.00 1.00
## embergerQ growingDegDays0 growingDegDays5 maxTempColdest
## Area 1.00 1.00 1.00 1.00
## minTempWarmest monthCountByTemp10 PETColdestQuarter
## Area 1.00 1.00 1.00
## PETDriestQuarter PETseasonality PETWarmestQuarter
## Area 1.00 1.00 1.00
## PETWettestQuarter thermicityIndex MaxD MinD
## Area 1.00 1.00 1.00 1.00
## [ reached getOption("max.print") -- omitted 37 rows ]
##
## To see confidence intervals of the correlations, print with the short=FALSE option
Διακρίνετε κάποιες μεταβλητές με τόσο μεγάλη συσχέτιση;
12.2.1 Αυτοματοποιημένος τρόπος
## ============================================================================##
## Correlation check
## ============================================================================##
unconc <- usdm::vifcor(ryukyus_new %>% as.data.frame() %>% dplyr::select(Area:MinD),
th = 0.7)
unconc <- unconc@results$Variables[unconc@results$VIF <= 5]
unconc <- droplevels(factor(unconc))
## ============================================================================##12.3 Οπτικοποίηση αποτελεσμάτων
Μπορούμε να οπτικοποιήσουμε τα αποτελέσματα μας και να διαπιστώσουμε εάν κάποια μεταβλητή δεν εμφανίζει κανονική κατανομή
## ============================================================================##
## Visualization
## ============================================================================##
car::scatterplotMatrix(ryukyus_new[, 5:10], spread = F, smoother.args = list(lty = 2),
main = "Scatter Plot Matrix")12.4 Κανονικοποίηση των δεδομένων μας
Καθώς τα δεδομένα μας δεν έχουν κανονική κατανομή, θα πρέπει να τα κανονικοποιήσουμε. Αυτό σημαίνει ότι θα πρέπει να τα λογαριθμήσουμε.
12.4.1 Λογαρίθμηση
Εφ’όσον έχουμε ελέγξει όλες τις μεταβλητές μας, μπορούμε να τις λογαριθμήσουμε4:
## ============================================================================##
## Log-transform the data
## ============================================================================##
ryukyus_log <- ryukyus_new %>% dplyr::select(Area, bio14, bio8, monthCountByTemp10,
MaxD, MinD) %>% mutate(SR = ryukyus_new$SR) %>% log10(.) %>% mutate(Island = ryukyus_new$Island,
PD = ryukyus_new$PD) %>% dplyr::select(Island, SR, PD, everything())
## ============================================================================##13 Κυρίως ανάλυση
13.1 Μοντελοποίηση - ISAR
Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 2 μοντέλα απλής γραμμικής παλινδρόμησης, ένα για την ταξινομική και ένα για την φυλογενετική ποικιλότητα σε σχέση με την έκταση:
## ============================================================================##
## ISAR models
## ============================================================================##
model1 <- lm(SR ~ Area, data = ryukyus_log)
model2 <- lm(PD ~ Area, data = ryukyus_log)
## ============================================================================##13.1.1 Περίληψη μοντέλων
Ας δούμε την περίληψη για κάθε ένα εκ των μοντέλων μας
## ============================================================================##
## Check the summary
## ============================================================================##
summary(model1)##
## Call:
## lm(formula = SR ~ Area, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.169622 -0.049189 -0.001934 0.072756 0.107531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.17927 0.08117 26.848 1.93e-13 ***
## Area 0.28263 0.03750 7.537 2.72e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08492 on 14 degrees of freedom
## Multiple R-squared: 0.8023, Adjusted R-squared: 0.7882
## F-statistic: 56.81 on 1 and 14 DF, p-value: 2.722e-06
## ============================================================================##
## Check the summary
## ============================================================================##
summary(model2)##
## Call:
## lm(formula = PD ~ Area, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02464 -0.50739 0.04005 0.60604 1.37789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.4217 0.6971 -3.474 0.00372 **
## Area 0.5013 0.3221 1.557 0.14188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7293 on 14 degrees of freedom
## Multiple R-squared: 0.1475, Adjusted R-squared: 0.08665
## F-statistic: 2.423 on 1 and 14 DF, p-value: 0.1419
ΠΡΟΣΟΧΗ
Μας ενδιαφέρει το adjusted R squared
13.1.2 Δημιουργία standardized residuals
Δημιουργούμε 2 αρχεία τα οποία περιέχουν τα standardized residuals των μοντέλων που δημιουργήσαμε στο προηγούμενο βήμα:
## ============================================================================##
## Create the residuals
## ============================================================================##
res_sr <- rstandard(model1)
res_pd <- rstandard(model2)
## ============================================================================##Μπορεί να γίνει σε ένα βήμα αντί για τρια:
## ============================================================================##
## Check the estimate and the p-value
## ============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% tidy() %>% full_join(., lm(PD ~
Area, data = ryukyus_log) %>% tidy()) %>% mutate(Metric = rep(c("SR", "PD"),
each = 2)) %>% dplyr::select(Metric, everything())
## ============================================================================##
## ============================================================================##
## Check the metrics
## ============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% glance() %>% full_join(., lm(PD ~
Area, data = ryukyus_log) %>% glance()) %>% mutate(Metric = c("SR", "PD")) %>%
dplyr::select(Metric, everything())
## ============================================================================##
## ============================================================================##
## Check the residuals
## ============================================================================##
resids <- lm(SR ~ Area, data = ryukyus_log) %>% augment() %>% mutate(Island = ryukyus_log$Island) %>%
rename(SR_res = .std.resid) %>% full_join(., lm(PD ~ Area, data = ryukyus_log) %>%
augment() %>% mutate(Island = ryukyus_log$Island) %>% rename(PD_res = .std.resid),
by = "Island") %>% dplyr::select(Island, SR, SR_res, PD, PD_res)
## ============================================================================##13.1.3 Στικτόγραμμα
Τώρα θα δημιουργήσουμε ένα στικτόγραμμα για να δούμε ποιό νησί αποτελεί θερμό σημείο φυτικής ποικιλότητας για τον συνολικό αριθμό των φυτικών taxa, καθώς και για τον αριθμό των ενδημικών taxa:
## ============================================================================##
## Plot
## ============================================================================##
ggplot(ryukyus_log, aes(res_sr, res_pd)) + geom_point() + geom_text(aes(label = Island),
fontface = "bold", size = 2.8, vjust = -1.2, hjust = 0.5) + geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) + theme(legend.position = "none") + xlab("SR") + ylab("PD")13.1.4 Αποθήκευση στικτογράμματος
Ας αποθηκεύσουμε την εικόνα αυτή:
13.2 Μοντελοποίηση - Πολλαπλή γραμμική παλινδρόμηση
13.2.1 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο enter
Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 2 μοντέλα πολλαπλής γραμμικής παλινδρόμησης, ένα για την ταξινομική και ένα για την φυλογενετική ποικιλότητα σε σχέση με όλες τις ανεξάρτητες, μη συγγραμικές μεταβλητές μας. Θα χρησιμοποιήσουμε αρχικά την μέθοδο enter:
## ============================================================================##
## The enter method
## ============================================================================##
enter1 <- lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)
enter2 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)
## ============================================================================##ΠΡΟΣΟΧΗ
Τα ανωτέρω μοντέλα δεν δείχνουν ποιές είναι οι μεταβλητές οι οποίες επηρεάζουν τα πρότυπα ποικιλότητας. Δείχνουν απλά τι ποσοστό εξηγείται εάν συμπεριλαμβάνονται όλες οι ανεξάρτητες μεταβλητές στο μοντέλο μας.
Ας δούμε τώρα την περίληψη για κάθε ένα εξ αυτών των μοντέλων:
## ============================================================================##
## Check the summary
## ============================================================================##
summary(enter1)##
## Call:
## lm(formula = SR ~ Area + bio14 + bio8 + monthCountByTemp10 +
## MaxD + MinD, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.135024 -0.049051 0.003223 0.031492 0.114373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.09862 2.01050 0.546 0.5980
## Area 0.27193 0.04050 6.715 8.7e-05 ***
## bio14 0.17172 0.48562 0.354 0.7318
## bio8 1.04033 0.59001 1.763 0.1117
## monthCountByTemp10 -1.52961 0.81899 -1.868 0.0946 .
## MaxD 0.02294 0.30523 0.075 0.9417
## MinD -0.09494 0.13400 -0.708 0.4966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08294 on 9 degrees of freedom
## Multiple R-squared: 0.8788, Adjusted R-squared: 0.7979
## F-statistic: 10.87 on 6 and 9 DF, p-value: 0.001092
##
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 +
## MaxD + MinD, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8957 -0.2648 -0.0292 0.4049 0.7872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -20.5226 14.1621 -1.449 0.1812
## Area 0.7156 0.2853 2.509 0.0334 *
## bio14 7.5665 3.4207 2.212 0.0543 .
## bio8 5.3023 4.1561 1.276 0.2340
## monthCountByTemp10 -5.9614 5.7691 -1.033 0.3284
## MaxD -0.3244 2.1501 -0.151 0.8834
## MinD -2.1573 0.9439 -2.285 0.0481 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5842 on 9 degrees of freedom
## Multiple R-squared: 0.6484, Adjusted R-squared: 0.4139
## F-statistic: 2.766 on 6 and 9 DF, p-value: 0.08269
13.2.2 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο της βηματικής πολλαπλής γραμμικής παλινδρόμησης
Με την μέθοδο αυτή θα μπορέσουμε να δούμε ποιές μεταβλητές εμπερικλείονται στο βέλτιστο μοντέλο για τον συνολικό αριθμό των ειδών
## ============================================================================##
## The stepAIC method
## ============================================================================##
stepm1 <- MASS::stepAIC(lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD +
MinD, data = ryukyus_log), scale = 0, direction = "both", trace = 1, keep = NULL,
steps = 1000, k = 2)
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD +
MinD, data = ryukyus_log), scale = 0, direction = "both", trace = 1, keep = NULL,
steps = 1000, k = 2)
summary(stepm1)Κάντε το ίδιο και για την φυλογενετική ποικιλότητα. Ποιές μεταβλητές εξηγούν και σε τι ποσοστό τα πρότυπα της φυλογενετικής ποικιλότητας;
ΠΡΟΣΟΧΗ: Μας ενδιαφέρει το adjusted R squared
## ============================================================================##
## The solution
## ============================================================================##
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD +
MinD, data = ryukyus_log), scale = 0, direction = "both", trace = 1, keep = NULL,
steps = 1000, k = 2)## Start: AIC=-12.41
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD
##
## Df Sum of Sq RSS AIC
## - MaxD 1 0.00777 3.0795 -14.3653
## - monthCountByTemp10 1 0.36444 3.4361 -12.6119
## <none> 3.0717 -12.4058
## - bio8 1 0.55550 3.6272 -11.7461
## - bio14 1 1.66988 4.7416 -7.4595
## - MinD 1 1.78274 4.8544 -7.0831
## - Area 1 2.14766 5.2194 -5.9234
##
## Step: AIC=-14.37
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD
##
## Df Sum of Sq RSS AIC
## - monthCountByTemp10 1 0.38000 3.4595 -14.5036
## <none> 3.0795 -14.3653
## - bio8 1 0.56038 3.6398 -13.6903
## + MaxD 1 0.00777 3.0717 -12.4058
## - bio14 1 1.73902 4.8185 -9.2021
## - MinD 1 1.85944 4.9389 -8.8071
## - Area 1 2.13989 5.2194 -7.9234
##
## Step: AIC=-14.5
## PD ~ Area + bio14 + bio8 + MinD
##
## Df Sum of Sq RSS AIC
## - bio8 1 0.33274 3.7922 -15.0342
## <none> 3.4595 -14.5036
## + monthCountByTemp10 1 0.38000 3.0795 -14.3653
## + MaxD 1 0.02333 3.4361 -12.6119
## - bio14 1 1.46668 4.9261 -10.8485
## - MinD 1 1.81014 5.2696 -9.7701
## - Area 1 2.38991 5.8494 -8.1001
##
## Step: AIC=-15.03
## PD ~ Area + bio14 + MinD
##
## Df Sum of Sq RSS AIC
## <none> 3.7922 -15.0342
## + bio8 1 0.33274 3.4595 -14.5036
## + monthCountByTemp10 1 0.15236 3.6398 -13.6903
## + MaxD 1 0.00319 3.7890 -13.0477
## - bio14 1 1.21835 5.0106 -12.5767
## - MinD 1 1.62428 5.4165 -11.3303
## - Area 1 2.38796 6.1802 -9.2199
##
## Call:
## lm(formula = PD ~ Area + bio14 + MinD, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99647 -0.28430 -0.05602 0.32450 0.81760
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -12.8387 6.5774 -1.952 0.0747 .
## Area 0.7480 0.2721 2.749 0.0176 *
## bio14 5.8728 2.9910 1.963 0.0732 .
## MinD -1.5102 0.6661 -2.267 0.0427 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5622 on 12 degrees of freedom
## Multiple R-squared: 0.5659, Adjusted R-squared: 0.4573
## F-statistic: 5.214 on 3 and 12 DF, p-value: 0.01554
## Start: AIC=-12.41
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD
##
## Df Sum of Sq RSS AIC
## - MaxD 1 0.00777 3.0795 -14.3653
## - monthCountByTemp10 1 0.36444 3.4361 -12.6119
## <none> 3.0717 -12.4058
## - bio8 1 0.55550 3.6272 -11.7461
## - bio14 1 1.66988 4.7416 -7.4595
## - MinD 1 1.78274 4.8544 -7.0831
## - Area 1 2.14766 5.2194 -5.9234
##
## Step: AIC=-14.37
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD
##
## Df Sum of Sq RSS AIC
## - monthCountByTemp10 1 0.38000 3.4595 -14.5036
## <none> 3.0795 -14.3653
## - bio8 1 0.56038 3.6398 -13.6903
## + MaxD 1 0.00777 3.0717 -12.4058
## - bio14 1 1.73902 4.8185 -9.2021
## - MinD 1 1.85944 4.9389 -8.8071
## - Area 1 2.13989 5.2194 -7.9234
##
## Step: AIC=-14.5
## PD ~ Area + bio14 + bio8 + MinD
##
## Df Sum of Sq RSS AIC
## - bio8 1 0.33274 3.7922 -15.0342
## <none> 3.4595 -14.5036
## + monthCountByTemp10 1 0.38000 3.0795 -14.3653
## + MaxD 1 0.02333 3.4361 -12.6119
## - bio14 1 1.46668 4.9261 -10.8485
## - MinD 1 1.81014 5.2696 -9.7701
## - Area 1 2.38991 5.8494 -8.1001
##
## Step: AIC=-15.03
## PD ~ Area + bio14 + MinD
##
## Df Sum of Sq RSS AIC
## <none> 3.7922 -15.0342
## + bio8 1 0.33274 3.4595 -14.5036
## + monthCountByTemp10 1 0.15236 3.6398 -13.6903
## + MaxD 1 0.00319 3.7890 -13.0477
## - bio14 1 1.21835 5.0106 -12.5767
## - MinD 1 1.62428 5.4165 -11.3303
## - Area 1 2.38796 6.1802 -9.2199
##
## Call:
## lm(formula = PD ~ Area + bio14 + MinD, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99647 -0.28430 -0.05602 0.32450 0.81760
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -12.8387 6.5774 -1.952 0.0747 .
## Area 0.7480 0.2721 2.749 0.0176 *
## bio14 5.8728 2.9910 1.963 0.0732 .
## MinD -1.5102 0.6661 -2.267 0.0427 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5622 on 12 degrees of freedom
## Multiple R-squared: 0.5659, Adjusted R-squared: 0.4573
## F-statistic: 5.214 on 3 and 12 DF, p-value: 0.01554
13.2.2.1 Ποιό μοντέλο είναι καλύτερο;
Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.
## ============================================================================##
## Compare the two models
## ============================================================================##
AIC(enter1, stepm1)Κάντε το ίδιο και για την φυλογενετική ποικιλότητα. Ποιό από τα δύο μοντέλα για κάθε κατηγορία ειδών είναι καλύτερο;
13.2.2.2 Δημιουργία μιας λειτουργίας για τον υπολογισμό της σχετικής σημαντικότητας των μεταβλητών
## ============================================================================##
## A function for the relative importance
## ============================================================================##
relweights <- function(fit, ...) {
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda^2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta^2)
rawwgt <- lambdasq %*% beta^2
import <- (rawwgt/rsquare) * 100
import <- as.data.frame(import)
row.names(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
import <- import[order(import), 1, drop = FALSE]
dotchart(import$Weights, labels = row.names(import), xlab = "% of R-Square",
pch = 19, main = "Relative Importance of Predictor Variables", sub = paste("Total R-Square=",
round(rsquare, digits = 3)), ...)
return(import)
}
## ============================================================================##13.2.2.3 Υπολογισμός της σχετικής σημαντικότητας των μεταβλητών
Μετά την περισπωμένη (~) βάζουμε όποιες μεταβλητές έχουν βγει σημαντικές από την βηματική πολλαπλή γραμμική παλινδρόμηση:
## ============================================================================##
## Check the relative importance
## ============================================================================##
rel.1 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)
relweights(rel.1, col = "blue")| Weights | |
|---|---|
| MaxD | 4.671829 |
| bio8 | 4.868050 |
| monthCountByTemp10 | 5.001657 |
| MinD | 22.456880 |
| bio14 | 31.196091 |
| Area | 31.805492 |
13.2.3 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο των ολικών υποσυνόλων
Μέχρι στιγμής, έχουμε δει δύο διαφορετικές μεθόδους, όσον αφορά την πολλαπλή γραμμική παλινδρόμηση. Όμως, δεν μπορούμε να είμαστε 100% σίγουροι ότι έχουμε φθάσει στο ορθό συμπέρασμα, καθότι ένα από τα μειονεκτήματα της βηματικής πολλαπλής παλινδρόμησης, είναι ότι μπορεί να “κολλήσει” σε ένα υποβέλτιστο πλατώ. Τον σκόπελο αυτό, μπορούμε να τον αποφύγουμε χρησιμοποιώντας την μέθοδο των ολικών υποσυνόλων:
## ============================================================================##
## The best subsets method
## ============================================================================##
bes.tot <- leaps::regsubsets(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD +
MinD, data = ryukyus_log, nbest = 4)
## ============================================================================##13.2.3.1 Απεικόνιση των αποτελεσμάτων του προηγούμενου βηματος
Το γράφημα αυτό θα μας δείξει ποιός είναι ο πραγματικά βέλτιστος συνδυασμός των ανεξάρτητων μεταβλητών:
## ============================================================================##
## Visualise
## ============================================================================##
plot(bes.tot, scale = "adjr2")## ============================================================================##
## Check the summary and the relative importance
## ============================================================================##
tot <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD, data = ryukyus_log)
summary(tot)##
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 +
## MinD, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.89706 -0.24475 -0.02862 0.37684 0.80405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -21.6286 11.5101 -1.879 0.0897 .
## Area 0.7130 0.2705 2.636 0.0249 *
## bio14 7.4248 3.1244 2.376 0.0389 *
## bio8 5.3227 3.9457 1.349 0.2071
## monthCountByTemp10 -5.6306 5.0687 -1.111 0.2926
## MinD -2.1181 0.8620 -2.457 0.0338 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5549 on 10 degrees of freedom
## Multiple R-squared: 0.6475, Adjusted R-squared: 0.4712
## F-statistic: 3.673 on 5 and 10 DF, p-value: 0.03798
| Weights | |
|---|---|
| bio8 | 5.224379 |
| monthCountByTemp10 | 5.706259 |
| MinD | 24.302432 |
| Area | 31.886016 |
| bio14 | 32.880914 |
13.2.4 Έλεγχος των αποτελεσμάτων της παλινδρόμησης
Τώρα θα ελέγξουμε κατ’αρχάς εάν οι μεταβλητές του βέλτιστου μοντέλου μας παραβιάζουν μια από τις αρχές της πολλαπλής γραμμικής παλινδρόμησης:
## Area bio14 bio8 monthCountByTemp10
## FALSE FALSE FALSE FALSE
## MinD
## FALSE
13.2.5 Ποιό μοντέλο είναι καλύτερο;
Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.
## ============================================================================##
## Compare the two models
## ============================================================================##
AIC(enter2, stepm2, tot)| df | AIC | |
|---|---|---|
| enter2 | 8 | 35.00027 |
| stepm2 | 5 | 32.37180 |
| tot | 7 | 33.04070 |
13.2.6 Έλεγχος πιθανοτήτων
Στο σημείο αυτό θα χρειαστεί να ελέγξουμε εάν η πιθανότητα κάθε μεταβλητής η οποία εμπερικλείεται στο τελικό μας μοντέλο είναι πραγματικά στατιστικώς σημαντική:
## ============================================================================##
## Let's have a look on the p-values for each of the IVs that have entered our
## final model
## ============================================================================##
summary(tot)##
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 +
## MinD, data = ryukyus_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.89706 -0.24475 -0.02862 0.37684 0.80405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -21.6286 11.5101 -1.879 0.0897 .
## Area 0.7130 0.2705 2.636 0.0249 *
## bio14 7.4248 3.1244 2.376 0.0389 *
## bio8 5.3227 3.9457 1.349 0.2071
## monthCountByTemp10 -5.6306 5.0687 -1.111 0.2926
## MinD -2.1181 0.8620 -2.457 0.0338 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5549 on 10 degrees of freedom
## Multiple R-squared: 0.6475, Adjusted R-squared: 0.4712
## F-statistic: 3.673 on 5 and 10 DF, p-value: 0.03798
## ============================================================================##
## Inside the parenthesis we include the above p-values
## ============================================================================##
p <- c(0.0249, 0.0389, 0.2071, 0.2926, 0.0338)
## ============================================================================##
## ============================================================================##
## Level of significance
## ============================================================================##
alp <- 0.05
## ============================================================================##
## ============================================================================##
## If the model's p-values are equal or smaller than the p-values we obtain from
## the following command, then they are ok
## ============================================================================##
p.adjust(p, method = "bonferroni", n = length(p))## [1] 0.1245 0.1945 1.0000 1.0000 0.1690
14 Θερμά σημεία βιοποικιλότητας
14.1 Εισαγωγή
Οι κλιμακούμενες απειλές για την βιοποικιλότητα, σε συνδυασμό με την υποχρηματοδότηση, έχουν ως αποτέλεσμα την εφαρμογή μιας διαδικασίας προτεραιοποίησης όσον αφορά τη λήψη αποφάσεων διαχειριστικού χαρακτήρα. Η διαδικασία αυτή επικεντρώνεται στην ιεράρχηση των ειδών, των πληθυσμών ή/και των οικοτόπων, με βάση τα οφέλη της βιοποικιλότητας, την δυνατότητα ανάκαμψης των ειδών και το κόστος επίτευξης του επιθυμητού στόχου (Farrier et al., 2007; Bottrill et al., 2008; Arponen, 2012; Jetz et al., 2014; Reece & Noss, 2014).
Ένα αποτελεσματικό εργαλείο για αυτόν τον σκοπό, είναι ο εντοπισμός θερμών σημείων ταξινομικής και φυλογενετικής βιοποικιλότητας.
Η κύρια στρατηγική για τη διατήρηση της βιοποικιλότητας έγκειται στον προσδιορισμό περιοχών για την προάσπιση των οργανισμών έναντι της εκάστοτε υφιστάμενης ή μελλοντικής απειλής (Margules & Pressey, 2000). Τα δίκτυα προστατευόμενων περιοχών (ΠΠ) αποτελούν τη ναυαρχίδα των πολιτικών διατήρησης της βιοποικιλότητας, ενώ διαδραματίζουν σημαίνοντα ρόλο στην προστασία της (Chape et al., 2005), καθώς η διατήρηση αρκετών ειδών είναι ισχυρώς εξαρτώμενη από τα δίκτυα αυτά (Watson et al., 2014). Παγκοσμίως, σχεδόν το 13% της χερσαίας έκτασης βρίσκεται υπό κάποιο καθεστώς προστασίας, ενώ το αντίστοιχο ποσοστό στις χώρες του ανεπτυγμένου κόσμου ανέρχεται στο ~12%, με εμφανή τάση αύξησης τα τελευταία χρόνια (Jenkins & Joppa, 2009). Για την αντιμετώπιση της παρατηρούμενης τάσης μείωσης της βιοποικιλότητας, η Συνθήκη για τη Βιολογική Ποικιλότητα, μέσω της Παγκόσμιας Στρατηγικής για τη Διατήρηση της Φυτικής και Ζωικής Ποικιλότητας, έχει θέσει ως στόχο την ανάσχεση της μείωσης αυτής (Jackson & Kennedy, 2009), δια της επίτευξης των Aichi Targets.
Η προτεραιοποίηση των περιοχών διατήρησης εδράζεται στον εντοπισμό θερμών σημείων βιοποικιλότητας (ΘΣΒ), τα οποία παραδοσιακά προσδιορίζονται βάσει ταξινομικών δεικτών [π.χ., αριθμός ειδών/ενδημικών, ποσοστό ενδημισμού, πιθανότητα εξαφάνισης – Forest et al. (2007)] και αποτελούν σημαντικά εργαλεία στην διαδικασία διαλογής όσον αφορά τη λήψη αποφάσεων διαχειριστικού χαρακτήρα. Ωστόσο, η προσέγγιση αυτή δεν λαμβάνει υπόψη ότι όλα τα είδη δεν είναι ισοδύναμα όσον αφορά την εξελικτική τους ιστορία ή/και τις λειτουργίες που επιτελούν (Cadotte & Jonathan Davies, 2010), ούτε το γεγονός ότι οι διαειδικές φυλογενετικές σχέσεις και κατ’ επέκταση η φυλογενετική ποικιλότητα της εκάστοτε περιοχής, αντικατοπτρίζουν τις οικολογικές, εξελικτικές και βιογεωγραφικές διαδικασίες με τις οποίες παράγεται, κατανέμεται και διατηρείται η βιοποικιλότητα (Cadotte & Jonathan Davies, 2010; Daru et al., 2019). Δεδομένου ότι ορισμένα είδη έχουν περισσότερο διακριτή εξελικτική ιστορία, η μη τυχαία εξαφάνιση ειδών μπορεί να οδηγήσει ορισμένους φυλογενετικούς κλάδους να απωλέσουν μεγαλύτερο ποσοστό ειδών σε σχέση με άλλους (Davies & Yessoufou, 2013). Καθώς είναι πλέον εφικτή η ενσωμάτωση φυλογενετικής πληροφορίας (Smith & Brown, 2018) στον προσδιορισμό και την προτεραιοποίηση των περιοχών διατήρησης, το γεγονός αυτό μπορεί να αποτελέσει ένα από τα σημαντικότερα βήματα προς την αποτελεσματικότερη διατήρηση και προστασία της βιοποικιλότητας (Tucker et al., 2019) σε διάφορες χωρικές και ταξινομικές κλίμακες (Pollock et al., 2017).
Στο φυλογενετικό επίπεδο, οι ισοδύναμοι των παραδοσιακών, ταξινομικών δεικτών για τον εντοπισμό ΘΣΒ, είναι η φυλογενετική ποικιλότητα [ΦΠ, sensu Faith (1992) – ισοδύναμη του αριθμού των ειδών], ο φυλογενετικός ενδημισμός [ΦΕ, sensu Rosauer et al. (2009) – ισοδύναμος του αριθμού των ενδημικών/ποσοστού ενδημισμού] και ο δείκτης EDGE [Evolutionary Distinct and Globally Endangered – EDGE, sensu Isaac et al. (2007) – ισοδύναμος της πιθανότητας εξαφάνισης]. Οι εν λόγω δείκτες ποσοτικοποιούν διαφορετικές πτυχές της εξελικτικής ποικιλότητας (Tucker et al., 2017).
Η ΦΠ είναι το άθροισμα του μήκους των κλάδων που συνδέουν ένα σύμπλεγμα ειδών με την ρίζα ενός φυλογενετικού δένδρου (Faith, 1992).
Ο ΦΕ ποσοτικοποιεί τον βαθμό στον οποίο ένα σημαντικό ποσοστό της ΦΠ εντοπίζεται αποκλειστικά εντός της εκάστοτε περιοχής μελέτης (Rosauer et al., 2009).
Ο δείκτης EDGE συνδυάζει την εξελικτική διακριτότητα (ED, η φυλογενετική απομόνωση ενός είδους) με το καθεστώς κινδύνου εξαφάνισης σε παγκόσμιο επίπεδο (GE) σύμφωνα με τα πρότυπα της IUCN (Isaac et al., 2007).
Οι δείκτες ΦΕ και EDGE αντιπροσωπεύουν σταθμισμένες παραλλαγές της ΦΠ σε γεωγραφική κλίμακα και καθεστώς απειλής, αντίστοιχα (Daru et al., 2019). Περιοχές με στατιστικώς σημαντικά υψηλό ή χαμηλό ΦΕ είναι κέντρα πάλαιο- ή νέο-ενδημισμού, αντίστοιχα (Mishler et al., 2014).
Ένας ακόμα σημαντικός δείκτης για τον ασφαλή εντοπισμό ΘΣΒ, είναι η σταθμισμένη ταξινομική ποικιλότητα (Corrected weighted endemism), καθώς δίνει περισσότερη έμφαση στα στενότοπα taxa (Crisp et al., 2001; Linder, 2001a, 2001b)
Συνεπώς, συνδυάζοντας τους ταξινομικούς και φυλογενετικούς αυτούς δείκτες, ο εντοπισμός των ΘΣΒ και κατ’ επέκταση η προτεραιοποίηση των ΠΠ αποκτά έναν πιο ολοκληρωμένο και ολιστικό χαρακτήρα που πιθανότατα επαρκεί για τη μακροπρόθεσμη διατήρηση της βιοποικιλότητας (Buerki et al., 2015), μειώνοντας το οικονομικό κόστος των μέτρων διατήρησης (Carta et al., 2019).
14.2 Εντοπισμός θερμών σημείων βιοποικιλότητας
Σε αυτό το σημείο θα χρησιμοποιήσουμε τα δεδομένα από την Ιταλία, τα οποία δημιουργήσαμε σε προηγούμενο εργαστήριο.
14.2.1 Εισαγωγή δεδομένων
Ας εισαγάγουμε τα δεδομένα για την Ιταλία.
Θυμηθείτε ότι θα χρειαστεί να αλλάξετε το file path και να χρησιμοποιήσετε τα δεδομένα από το Github.
## ============================================================================##
## Load the biotic cleaned data, as well as the spatial data for Italy
## ============================================================================##
data_it_cleaned <- readRDS("Italian data cleaned.rds")
Italy <- read_sf("Italin provinces.shp") %>% st_transform(4326) %>% as_Spatial()14.2.2 Προκαταρτικά βήματα
Θα χρειαστεί να καθαρίσουμε τα ονόματα των taxa5 και στη συνέχεια να κρατήσουμε μόνο τα έγκυρα ονόματα.
## ============================================================================##
## First load the excel file with the correct species names
## ============================================================================##
nms <- readxl::read_excel("italian_names.xlsx")
## ============================================================================##
## ============================================================================##
## Keep only one instance per taxon from our occurrence data
## ============================================================================##
nms_flags <- data_it_cleaned %>% distinct(scientificName) %>% mutate(Taxon = nms$scientificName,
action = nms$action) %>% as_tibble()
## ============================================================================##Ελέγξτε εάν υπάρχει διαφορά μεταξύ των δύο αρχείων.
14.2.3 Καθαρισμός δεδομένων
Στη συνέχεια θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει τις συντεταγμένες για τα σωστά taxa:
## ============================================================================##
## First create the coords for the correct taxa
## ============================================================================##
coords <- data_it_cleaned %>% full_join(nms_flags) %>% as_tibble() %>% filter(!str_detect(action,
"DEL")) %>% dplyr::select(decimalLongitude, decimalLatitude, Taxon)
## ============================================================================##14.2.4 Δημιουργία μήτρας παρουσίας/απουσίας
Και ακολούθως ένα αρχείο το οποίο θα αποτελέσει την βάση των αναλύσεων μας:
## ============================================================================##
## Then create the sparse matrix
## ============================================================================##
pt <- points2comm(dat = coords, mask = Greece, res = 0.04166667, lon = "decimalLongitude",
lat = "decimalLatitude", species = "Taxon")Ας δούμε πως είναι:
## 5 x 5 sparse Matrix of class "dgCMatrix"
## Abies_alba Abies_nebrodensis Abies_nordmanniana Abutilon_theophrasti
## v10000 . . . .
## v10003 . . . .
## v10004 . . . .
## v10007 . . . .
## v10008 . . . .
## Acacia_dealbata
## v10000 .
## v10003 .
## v10004 .
## v10007 .
## v10008 .
Και ας το αποθηκεύσουμε:
14.2.5 Κατασκευή φυλογενετικού δένδρου
Ας κατασκευάσουμε το φυλογενετικό δένδρο για τα taxa τα οποία απαντώνται στην Ιταλία.
Υπάρχει πιθανότητα αναλόγως των δυνατοτήτων του Η/Υ σας, η εντολή αυτή να διαρκέσει μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτής με την εντολή:
italy_tree <- readRDS(Italian phylogenetic tree.rds.
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
italy_tree <- data_it_cleaned %>%
full_join(nms_flags) %>%
as_tibble() %>%
filter(!str_detect(action, 'DEL')) %>%
dplyr::select(family, genus, Taxon) %>%
dplyr::rename(species = Taxon) %>%
distinct(species, .keep_all = T) %>%
get_tree(sp_list = .,
tree = tree_plant_otl, # either
taxon = "plant", # or
scenario = "S1",
show_grafted = FALSE)
##============================================================================##
## Save it to disk
##============================================================================##
saveRDS(italy_tree, 'Italian phylogenetic tree.rds')
##============================================================================##14.2.6 Κυρίως ανάλυση
Τώρα μπορούμε να δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει πληροφορία σχετικά με το ποια κελιά αποτελούν θερμά σημεία σταθμισμένης και μη-σταθμισμένης ταξινομικής και φυλογενετικής ποικιλότητας:
##============================================================================##
## Locate the hotspots
##============================================================================##
hot_cold_spots <- pt$poly_shp %>%
## Convert to sf
st_as_sf() %>%
## Create the hot- and cold-spots
mutate(Cold_SR = coldspots(richness,
prob = 5),
Hot_SR = hotspots(richness,
prob = 5),
## Corrected weighted endemism
Cold_CWE = coldspots(weighted_endemism(pt$comm_dat),
prob = 5),
Hot_CWE = hotspots(weighted_endemism(pt$comm_dat),
prob = 5),
## Phylogenetic diversity
Cold_PD = coldspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
Hot_PD = hotspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
## Phylogenetic endemism
Cold_PE = coldspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
Hot_PE = hotspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5)) %>%
mutate(all_hot = ifelse(Hot_SR == 1 & Hot_CWE == 1 & Hot_PD == 1 & Hot_PE == 1, 'HOT',
ifelse(Cold_SR == 1 & Cold_CWE == 1 & Cold_PD == 1 & Cold_PE == 1, 'COLD', 'Neither'))) %>%
mutate(cwe_sr_hot = ifelse(Hot_SR == 1 & Hot_CWE == 1, 'HOT',
ifelse(Cold_SR == 1 & Cold_CWE == 1, 'COLD', 'Neither')))
##============================================================================##14.2.7 Οπτικοποίηση
Ας δημιουργήσουμε ένα νέο αρχείο το οποίο θα χρησιμοποιήσουμε για την οπτικοποίηση των αποτελέσματων μας και το οποίο θα περιέχει μια νέα μεταβλητή, η οποία θα υποδεικνύει ποια κελιά αποτελούν θερμά ή ψυχρά σημεία βιοποικιλότητας για όλους τους δείκτες που χρησιμοποιήσαμε:
## ============================================================================##
## Create a new variable
## ============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>% as_Spatial()
## ============================================================================##
## ============================================================================##
## Quick plot
## ============================================================================##
plot(sp_hot_cold_spots, border = "grey", col = "lightgrey", main = "Diversity Hotspots and Coldspots")
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == "HOT"), ], col = "red",
add = TRUE, border = NA)
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == "COLD"), ], col = "blue",
add = TRUE, border = NA)
legend("bottomleft", fill = c("red", "blue"), legend = c("hotspots", "coldspots"),
bty = "n", inset = 0.092)## ============================================================================##
## Second quick plot
## ============================================================================##
hot_cold_spots %>% ggplot() + geom_sf(aes(fill = cwe_sr_hot)) + labs(title = "Diversity hotspots in Italy",
subtitle = "Source: GBIF") + scale_fill_manual(values = c("steelblue1", "red",
"gray90"), name = "Diversity hotspots", labels = c("COLD", "HOT", "Neither"))Απεικονίστε τα θερμά σημεία βιοποικιλότητας σύμφωνα με τους δείκτες PD & PE.
15 Appendix: R code
##===========================================================================##
## Do NOT run this code
##===========================================================================##
## Global options
options(Encoding="UTF-8")
library(knitr)
library(rmdformats)
## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
cache=TRUE,
prompt=FALSE,
tidy=TRUE,
# comment=NA,
message=FALSE,
warning=FALSE,
eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
htmltools::img(src = knitr::image_uri(file.path("G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/HBS Seminar/Logo_HBS/ebe-LOGO-005-300_1_70.jpg")),
alt = 'logo',
style = 'position:absolute; top:0; right:0; padding:5px;')
##===========================================================================##
## Load the libraries
##===========================================================================##
library(tidyverse)
library(raster)
library(sf)
library(ape)
library(phyloregion)
library(picante)
library(rtrees)
library(ggtree)
library(exactextractr)
library(caret)
library(psych)
library(PhyloMeasures)
library(phyloregion)
library(Matrix)
library(broom)
##===========================================================================##
##===========================================================================##
## Load our data
##===========================================================================##
plants <- read.csv("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.csv", h = T, sep = ";", row.names = 1)
##===========================================================================##
##===========================================================================##
## Check their dimensions
##===========================================================================##
dim(plants)
##===========================================================================##
## Check their dimensions
##===========================================================================##
str(plants)
pre {
max-height: 400px;
overflow-y: auto;
}
pre[class] {
max-height: 400px;
}
##===========================================================================##
## Convert the P/A matrix to SR
##===========================================================================##
ryukyus <- rowSums(plants) %>%
enframe() %>%
dplyr::rename(Island = name,
SR = value)
ryukyus
##===========================================================================##
## Create a vector containing just the file paths
##===========================================================================##
shps <- fs::dir_ls("C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/",
glob = "*.shp")
##===========================================================================##
##===========================================================================##
## Create a vector containing just the island names
##===========================================================================##
island_nms <- shps %>%
str_replace_all('C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/',
'') %>%
str_replace_all('\\.shp', '')
##===========================================================================##
##===========================================================================##
## Bulk load the shapefiles
##===========================================================================##
ryukyus_shps <- shps %>%
purrr::map(read_sf) %>%
lapply(., dplyr::select, NAME_1, geometry) %>%
do.call(rbind, .) %>%
mutate(Island = island_nms) %>%
dplyr::select(Island, geometry)
##===========================================================================##
##===========================================================================##
## Load, crop and mask the rasters
##===========================================================================##
preds <- stack(stack(list.files(path = "I:/SDM DATA/WORLDCLIM/Current/bio_2-5m_bil/",
pattern = '.bil',
recursive = T,
full.names = T)) %>%
crop(., ryukyus_shps, snap = 'in') %>%
mask(., ryukyus_shps),
stack(list.files(path = "I:/SDM DATA/Envirem_Predictors/Eurasia_current_2.5arcmin_geotiff/",
pattern = '.tif',
recursive = T,
full.names = T)) %>%
crop(., ryukyus_shps, snap = 'in') %>%
mask(., ryukyus_shps))
##===========================================================================##
##===========================================================================##
## Load, crop and mask the altitudinal raster
##===========================================================================##
altitude <- raster("I:/SDM DATA/Altitude/elevation_1KMmd_GMTEDmd.tif") %>%
crop(., ryukyus_shps, snap = 'in') %>%
mask(., ryukyus_shps)
##===========================================================================##
##===========================================================================##
## Extract abiotic info
##===========================================================================##
preds_ryukyus <- exact_extract(preds, ryukyus_shps, 'median') %>%
mutate(Island = ryukyus_shps$Island) %>%
full_join(.,
exact_extract(altitude, ryukyus_shps, 'median') %>%
enframe() %>%
mutate(Island = ryukyus_shps$Island) %>%
rename(Elevation = value) %>%
dplyr::select(Island, Elevation)) %>%
dplyr::select(Island, everything()) %>%
as_tibble()
##============================================================================##
## Give them proper names
##============================================================================##
col_names <- preds_ryukyus %>%
colnames() %>%
enframe() %>%
mutate(value = str_replace_all(value, 'median\\.', '')) %>%
mutate(value = str_replace_all(value, 'current_2\\.5arcmin_', '')) %>%
pull(value)
##===========================================================================##
##============================================================================##
## Overwrite names and remove Row 1
##============================================================================##
vars <- preds_ryukyus %>%
set_names(col_names) %>%
dplyr::select(Island, everything()) %>%
as_tibble()
##============================================================================##
## Transform the missing data
##============================================================================##
imputed_vars <- predict(preProcess(as.data.frame(vars[, 2:36]),
method = c('medianImpute')),
vars)
##============================================================================##
## First drop any rows with NA values
##============================================================================##
drop_vars <- vars %>% drop_na()
##============================================================================##
##============================================================================##
## Then check which islands have NA values in any of the variables
##============================================================================##
setdiff(imputed_vars$Island, drop_vars$Island)
##============================================================================##
## Create a symbol to filter out all the Islands having NA values
##============================================================================##
'%!in%' <- function(x,y)!('%in%'(x,y))
##============================================================================##
##============================================================================##
## Now create a temporary presence/absence matrix removing any problematic islands
##============================================================================##
plants_new <- plants %>%
mutate(Island = rownames(.)) %>%
as_tibble() %>%
dplyr::select(Island, everything()) %>%
pivot_longer(!Island, names_to = 'Taxon',
values_to = 'PA') %>%
filter(Island %!in% setdiff(imputed_vars$Island,
drop_vars$Island)) %>%
pivot_wider(names_from = Taxon,
values_from = PA)
##============================================================================##
##============================================================================##
## See which taxa are left after removing the problematic islands
##============================================================================##
keepers <- colSums(plants_new[, -1]) %>%
enframe() %>%
filter(value > 0)
##============================================================================##
##============================================================================##
## Finally create the new PA matrix
##============================================================================##
plants_new_pa <- plants_new %>%
pivot_longer(!Island, names_to = 'Taxon',
values_to = 'PA') %>%
filter(Taxon %in% keepers$name) %>%
pivot_wider(names_from = Taxon,
values_from = PA)
##============================================================================##
datatable(plants, options = list(pagelength = 10))
jpn_tree <- readRDS('G:/Academic/Lessons/EMB/Labs/Labs/16 Ryukyus tree.rds')
##============================================================================##
## Plot the tree
##============================================================================##
groupInfo <- split(jpn_tree$tip.label,
gsub("_\\w+", "",
jpn_tree$tip.label))
groups_jpn <- groupOTU(jpn_tree,
groupInfo)
ggtree(groups_jpn, aes(color = group),
layout = 'circular') +
geom_tiplab2(size = 0.5, align = T)
##============================================================================##
## Combine the trees with the P/A data
##============================================================================##
#check for mismatches/missing species
combined_japan <- picante::match.phylo.comm(jpn_tree,plants_new_pa[, 2:1800])
##============================================================================##
## The resulting object is a list with $phy and $comm elements
## replace our original data with the sorted/matched data
##============================================================================##
phy_japan <- combined_japan$phy
comm_japan <- combined_japan$comm
##============================================================================##
##============================================================================##
## Compute phylogenetic alpha diversity
##============================================================================##
## For the estimation of SES values
japan_pd <- pd.query(phy_japan, comm_japan, TRUE)
## For the p-values
japan_pd_p <- pd.pvalues(phy_japan, comm_japan, reps = 1000)
##============================================================================##
##============================================================================##
## Create the final dataset
##============================================================================##
ryukyus_new <- ryukyus %>%
filter(Island %in% plants_new_pa$Island) %>%
full_join(., drop_vars) %>%
full_join(., ryukyus_shps) %>%
filter(Island %!in% setdiff(imputed_vars$Island,
drop_vars$Island)) %>%
st_as_sf() %>%
## Lambert Azimuthal Equal Area projection
st_transform('+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs') %>%
## Create the Area and PD variables
mutate(Area = as.numeric(units::set_units(st_area(.), km^2)),
PD = japan_pd,
PD_P = japan_pd_p) %>%
## Transform back to WGS84
st_transform(4326) %>%
dplyr::select(Island, SR, PD, PD_P, Area, Elevation, everything()) %>%
## It has zero-variance
dplyr::select(-aridityIndexThornthwaite) %>%
## Drop the geometry
st_drop_geometry() %>%
## Join these data with the Maximum and Minimum distance between each island
## pair
inner_join(.,
## Calculate the Maximum distance between each island pair
apply((st_distance(ryukyus_shps %>%
st_centroid(),
ryukyus_shps %>%
st_centroid())/1e+3) %>%
as_tibble() %>%
set_names(ryukyus_shps$Island) %>%
units::drop_units() %>%
na_if(., 0), 1, max, na.rm = T) %>%
enframe() %>%
rename(MaxD = value) %>%
mutate(Island = ryukyus_shps$Island) %>%
dplyr::select(Island, MaxD) %>%
full_join(.,
## Calculate the Minimum distance between each island
## pair
apply((st_distance(ryukyus_shps %>%
st_centroid(),
ryukyus_shps %>%
st_centroid())/1e+3) %>%
as_tibble() %>%
set_names(ryukyus_shps$Island) %>%
units::drop_units() %>%
na_if(., 0), 1, min, na.rm = T) %>%
enframe() %>%
rename(MinD = value) %>%
mutate(Island = ryukyus_shps$Island) %>%
dplyr::select(Island, MinD)))
##============================================================================##
ryukyus_new
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
library(rtrees)
jpn_tree <- keepers %>%
separate(name, into = c('genus', 'binomial'), sep = '_', remove = F) %>%
rename(species = name) %>%
dplyr::select(species, genus) %>%
get_tree(sp_list = .,
tree = tree_plant_otl, # either
taxon = "plant", # or
scenario = "S1",
show_grafted = FALSE)
##============================================================================##
##============================================================================##
## Normality check
##============================================================================##
normality <- lapply(ryukyus_new[,2:42], shapiro.test)
##============================================================================##
normality
##============================================================================##
## Correlation check
##============================================================================##
corr.test(ryukyus_new[,5:42], method = "spearman")
##============================================================================##
## Visualization
##============================================================================##
car::scatterplotMatrix(ryukyus_new[,5:10],
spread = F,
smoother.args = list(lty = 2),
main = "Scatter Plot Matrix")
##============================================================================##
## Correlation check
##============================================================================##
unconc <- usdm::vifcor(ryukyus_new %>%
as.data.frame() %>%
dplyr::select(Area:MinD),
th = 0.7)
unconc <- unconc@results$Variables[unconc@results$VIF <= 5]
unconc <- droplevels(factor(unconc))
##============================================================================##
##============================================================================##
## Log-transform the data
##============================================================================##
ryukyus_log <- ryukyus_new %>%
dplyr::select(Area, bio14, bio8,
monthCountByTemp10,
MaxD, MinD) %>%
mutate(SR = ryukyus_new$SR) %>%
log10(.) %>%
mutate(Island = ryukyus_new$Island,
PD = ryukyus_new$PD) %>%
dplyr::select(Island, SR, PD, everything())
##============================================================================##
##============================================================================##
## ISAR models
##============================================================================##
model1 <- lm(SR ~ Area, data = ryukyus_log)
model2 <- lm(PD ~ Area, data = ryukyus_log)
##============================================================================##
##============================================================================##
## Check the summary
##============================================================================##
summary(model1)
##============================================================================##
## Check the summary
##============================================================================##
summary(model2)
##============================================================================##
## Create the residuals
##============================================================================##
res_sr <- rstandard(model1)
res_pd <- rstandard(model2)
##============================================================================##
##============================================================================##
## Check the estimate and the p-value
##============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% tidy() %>%
full_join(., lm(PD ~ Area, data = ryukyus_log) %>% tidy()) %>%
mutate(Metric = rep(c('SR','PD'), each = 2)) %>%
dplyr::select(Metric, everything())
##============================================================================##
##============================================================================##
## Check the metrics
##============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% glance() %>%
full_join(., lm(PD ~ Area, data = ryukyus_log) %>% glance()) %>%
mutate(Metric = c('SR','PD')) %>%
dplyr::select(Metric, everything())
##============================================================================##
##============================================================================##
## Check the residuals
##============================================================================##
resids <- lm(SR ~ Area,
data = ryukyus_log) %>%
augment() %>%
mutate(Island = ryukyus_log$Island) %>%
rename(SR_res = .std.resid) %>%
full_join(., lm(PD ~ Area, data = ryukyus_log) %>%
augment() %>%
mutate(Island = ryukyus_log$Island) %>%
rename(PD_res = .std.resid),
by = 'Island') %>%
dplyr::select(Island, SR, SR_res, PD, PD_res)
##============================================================================##
##============================================================================##
## Plot
##============================================================================##
ggplot(ryukyus_log, aes(res_sr, res_pd)) +
geom_point() +
geom_text(aes(label = Island),
fontface = "bold",
size = 2.8,
vjust = -1.2,
hjust = 0.5) +
geom_hline(yintercept = 0.0) +
geom_vline(xintercept = 0.0) +
theme(legend.position = "none") +
xlab("SR") +
ylab("PD")
##============================================================================##
## Quick save
##============================================================================##
ggsave("SR-PD.png", width = 20, height = 20, units = "cm", dpi = 600)
##============================================================================##
## The enter method
##============================================================================##
enter1 <- lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)
enter2 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)
##============================================================================##
##============================================================================##
## Check the summary
##============================================================================##
summary(enter1)
summary(enter2)
##============================================================================##
## The solution
##============================================================================##
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
scale = 0,
direction = "both",
trace = 1,
keep = NULL,
steps = 1000,
k = 2)
summary(stepm2)
##============================================================================##
## The stepAIC method
##============================================================================##
stepm1 <- MASS::stepAIC(lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
scale = 0,
direction = "both",
trace = 1,
keep = NULL,
steps = 1000,
k = 2)
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
scale = 0,
direction = "both",
trace = 1,
keep = NULL,
steps = 1000,
k = 2)
summary(stepm1)
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
scale = 0,
direction = "both",
trace = 1,
keep = NULL,
steps = 1000,
k = 2)
summary(stepm2)
##============================================================================##
## Compare the two models
##============================================================================##
AIC(enter1, stepm1)
##============================================================================##
## A function for the relative importance
##============================================================================##
relweights <- function(fit,...){
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda ^ 2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta ^ 2)
rawwgt <- lambdasq %*% beta ^ 2
import <- (rawwgt / rsquare) * 100
import <- as.data.frame(import)
row.names(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
import <- import[order(import),1, drop=FALSE]
dotchart(import$Weights, labels=row.names(import),
xlab="% of R-Square", pch=19,
main="Relative Importance of Predictor Variables",
sub=paste("Total R-Square=", round(rsquare, digits=3)),
...)
return(import)
}
##============================================================================##
##============================================================================##
## Check the relative importance
##============================================================================##
rel.1 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)
relweights(rel.1, col="blue")
##============================================================================##
## The best subsets method
##============================================================================##
bes.tot <- leaps::regsubsets(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD,
data = ryukyus_log,
nbest = 4)
##============================================================================##
##============================================================================##
## Visualise
##============================================================================##
plot(bes.tot, scale = "adjr2")
##============================================================================##
## Check the summary and the relative importance
##============================================================================##
tot <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD, data = ryukyus_log)
summary(tot)
relweights(tot, col="blue")
sqrt(car::vif(tot)) > 2
##============================================================================##
## Compare the two models
##============================================================================##
AIC(enter2, stepm2, tot)
##============================================================================##
## Load the biotic cleaned data, as well as the spatial data for Italy
##============================================================================##
data_it_cleaned <- readRDS('Italian data cleaned.rds')
Italy <- read_sf('Italin provinces.shp') %>%
st_transform(4326) %>%
as_Spatial()
##============================================================================##
## First load the excel file with the correct species names
##============================================================================##
nms <- readxl::read_excel('italian_names.xlsx')
##============================================================================##
##============================================================================##
## Keep only one instance per taxon from our occurrence data
##============================================================================##
nms_flags <- data_it_cleaned %>%
distinct(scientificName) %>%
mutate(Taxon = nms$scientificName,
action = nms$action) %>%
as_tibble()
##============================================================================##
##============================================================================##
## First create the coords for the correct taxa
##============================================================================##
coords <- data_it_cleaned %>%
full_join(nms_flags) %>%
as_tibble() %>%
filter(!str_detect(action, 'DEL')) %>%
dplyr::select(decimalLongitude,
decimalLatitude,
Taxon)
##============================================================================##
pt <- readRDS('Sparse matrix and relevant polygons for Italy.rds')
##============================================================================##
## Then create the sparse matrix
##============================================================================##
pt <- points2comm(dat = coords,
mask = Greece,
res = 0.04166667,
lon = "decimalLongitude",
lat = "decimalLatitude",
species = "Taxon")
head(pt[[1]][1:5, 1:5])
##============================================================================##
## Save it
##============================================================================##
saveRDS(pt, 'Sparse matrix and relevant polygons for Italy.rds')
italy_tree <- readRDS('Italian phylogenetic tree.rds')
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
italy_tree <- data_it_cleaned %>%
full_join(nms_flags) %>%
as_tibble() %>%
filter(!str_detect(action, 'DEL')) %>%
dplyr::select(family, genus, Taxon) %>%
dplyr::rename(species = Taxon) %>%
distinct(species, .keep_all = T) %>%
get_tree(sp_list = .,
tree = tree_plant_otl, # either
taxon = "plant", # or
scenario = "S1",
show_grafted = FALSE)
##============================================================================##
## Save it to disk
##============================================================================##
saveRDS(italy_tree, 'Italian phylogenetic tree.rds')
##============================================================================##
##============================================================================##
## Let's have a look on the p-values for each of the IVs that have entered
## our final model
##============================================================================##
summary(tot)
##============================================================================##
## Inside the parenthesis we include the above p-values
##============================================================================##
p <- c(0.0249, 0.0389, 0.2071, 0.2926, 0.0338)
##============================================================================##
##============================================================================##
## Level of significance
##============================================================================##
alp <- 0.05
##============================================================================##
##============================================================================##
## If the model's p-values are equal or smaller than the p-values we obtain
## from the following command, then they are ok
##============================================================================##
p.adjust(p, method = "bonferroni", n = length(p))
##============================================================================##
## Create a new variable
##============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>%
as_Spatial()
##============================================================================##
##============================================================================##
## Quick plot
##============================================================================##
plot(sp_hot_cold_spots,
border = "grey",
col = "lightgrey",
main = "Diversity Hotspots and Coldspots")
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == 'HOT'), ],
col = "red",
add = TRUE,
border = NA)
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == 'COLD'), ],
col = "blue",
add = TRUE,
border = NA)
legend("bottomleft",
fill = c("red", "blue"),
legend = c("hotspots", "coldspots"),
bty = "n",
inset = .092)
##============================================================================##
## Second quick plot
##============================================================================##
hot_cold_spots %>%
ggplot() +
geom_sf(aes(fill = cwe_sr_hot )) +
labs(title = 'Diversity hotspots in Italy',
subtitle = 'Source: GBIF') +
scale_fill_manual(values = c('steelblue1','red', 'gray90'),
name = "Diversity hotspots",
labels = c("COLD", "HOT", "Neither"))
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##Βιβλιογραφία
Arponen, A. (2012) Prioritizing species for conservation planning. Biodiversity and Conservation, 21, 875–893.
Bottrill, M.C., Joseph, L.N., Carwardine, J., Bode, M., Cook, C., Game, E.T., Grantham, H., Kark, S., Linke, S., McDonald-Madden, E., & others (2008) Is conservation triage just smart decision making? Trends in ecology & evolution, 23, 649–654.
Buerki, S., Callmander, M.W., Bachman, S., Moat, J., Labat, J.-N., & Forest, F. (2015) Incorporating evolutionary history into conservation planning in biodiversity hotspots. Philosophical Transactions of the Royal Society B: Biological Sciences, 370, 20140014.
Cadotte, M.W. & Jonathan Davies, T. (2010) Rarest of the rare: Advances in combining evolutionary distinctiveness and scarcity to inform conservation at biogeographical scales. Diversity and Distributions, 16, 376–385.
Carta, A., Gargano, D., Rossi, G., Bacchetta, G., Fenu, G., Montagnani, C., Abeli, T., Peruzzi, L., & Orsenigo, S. (2019) Phylogenetically informed spatial planning as a tool to prioritise areas for threatened plant conservation within a mediterranean biodiversity hotspot. Science of The Total Environment, 665, 1046–1052.
Chape, S., Harrison, J., Spalding, M., & Lysenko, I. (2005) Measuring the extent and effectiveness of protected areas as an indicator for meeting global biodiversity targets. Philosophical Transactions of the Royal Society B: Biological Sciences, 360, 443–455.
Crisp, M.D., Laffan, S., Linder, H.P., & Monro, A. (2001) Endemism in the australian flora. Journal of Biogeography, 28, 183–198.
Daru, B.H., Roux, P.C. le, Gopalraj, J., Park, D.S., Holt, B.G., & Greve, M. (2019) Spatial overlaps between the global protected areas network and terrestrial hotspots of evolutionary diversity. Global Ecology and Biogeography, 28, 757–766.
Davies, T.J. & Yessoufou, K. (2013) Revisiting the impacts of non-random extinction on the tree-of-life. Biology letters, 9, 20130343.
Faith, D.P. (1992) Conservation evaluation and phylogenetic diversity. Biological conservation, 61, 1–10.
Farrier, D., Whelan, R., & Mooney, C. (2007) Threatened species listing as a trigger for conservation action. Environmental Science & Policy, 10, 219–229.
Forest, F., Grenyer, R., Rouget, M., Davies, T.J., Cowling, R.M., Faith, D.P., Balmford, A., Manning, J.C., Procheş, Ş., Bank, M. van der, & others (2007) Preserving the evolutionary potential of floras in biodiversity hotspots. Nature, 445, 757–760.
Good, R. (1845) The geography of the flowering plants. Longsmans,
Isaac, N.J., Turvey, S.T., Collen, B., Waterman, C., & Baillie, J.E. (2007) Mammals on the edge: Conservation priorities based on threat and phylogeny. PloS one, 2, e296.
Jackson, P.W. & Kennedy, K. (2009) The global strategy for plant conservation: A challenge and opportunity for the international community. Trends in plant science, 14, 578–580.
Jenkins, C.N. & Joppa, L. (2009) Expansion of the global terrestrial protected area system. Biological conservation, 142, 2166–2174.
Jetz, W., Thomas, G.H., Joy, J.B., Redding, D.W., Hartmann, K., & Mooers, A.O. (2014) Global distribution and conservation of evolutionary distinctness in birds. Current biology, 24, 919–930.
Kira, T. (1991) Forest ecosystems of east and southeast asia in a global perspective. Ecological Research, 6, 185–200.
Li, D., Trotta, L., Marx, H.E., Allen, J.M., Sun, M., Soltis, D.E., Soltis, P.S., Guralnick, R.P., & Baiser, B. (2019) For common community phylogenetic analyses, go ahead and use synthesis phylogenies. Ecology, 100, e02788.
Linder, H. (2001a) On areas of endemism, with an example from the african restionaceae. Systematic biology, 50, 892–912.
Linder, H. (2001b) Plant diversity and endemism in sub-saharan tropical africa. Journal of Biogeography, 28, 169–182.
Maekawa, F. (1974) Origin and characteristics of japan’s flora. The Flora and Vegetation of Japan., 33–86.
Margules, C.R. & Pressey, R.L. (2000) Systematic conservation planning. Nature, 405, 243–253.
Mishler, B.D., Knerr, N., González-Orozco, C.E., Thornhill, A.H., Laffan, S.W., & Miller, J.T. (2014) Phylogenetic measures of biodiversity and neo-and paleo-endemism in australian acacia. Nature Communications, 5, 1–10.
Nakamura, K., Suwa, R., Denda, T., & Yokota, M. (2009) Geohistorical and current environmental influences on floristic differentiation in the ryukyu archipelago, japan. Journal of Biogeography, 36, 919–928.
Ota, H. (1998) Geographic patterns of endemism and speciation in amphibians and reptiles of the ryukyu archipelago, japan, with special reference to their paleogeographical implications. Researches on Population Ecology, 40, 189–204.
Pollock, L.J., Thuiller, W., & Jetz, W. (2017) Large conservation gains possible for global biodiversity facets. Nature, 546, 141–144.
Reece, J.S. & Noss, R.F. (2014) Prioritizing species by conservation value and vulnerability: A new index applied to species threatened by sea-level rise and other risks in florida. Natural Areas Journal, 34, 31–45.
Rosauer, D., Laffan, S.W., Crisp, M.D., Donnellan, S.C., & Cook, L.G. (2009) Phylogenetic endemism: A new approach for identifying geographical concentrations of evolutionary history. Molecular ecology, 18, 4061–4072.
Smith, S.A. & Brown, J.W. (2018) Constructing a broadly inclusive seed plant phylogeny. American journal of botany, 105, 302–314.
Triantis, K.A., Guilhaumon, F., & Whittaker, R.J. (2012) The island species–area relationship: Biology and statistics. Journal of Biogeography, 39, 215–231.
Tucker, C.M., Aze, T., Cadotte, M.W., Cantalapiedra, J.L., Chisholm, C., Dı́az, S., Grenyer, R., Huang, D., Mazel, F., Pearse, W.D., & others (2019) Assessing the utility of conserving evolutionary history. Biological Reviews, 94, 1740–1760.
Tucker, C.M., Cadotte, M.W., Carvalho, S.B., Davies, T.J., Ferrier, S., Fritz, S.A., Grenyer, R., Helmus, M.R., Jin, L.S., Mooers, A.O., & others (2017) A guide to phylogenetic metrics for conservation, community ecology and macroecology. Biological Reviews, 92, 698–715.
Watson, J.E., Dudley, N., Segan, D.B., & Hockings, M. (2014) The performance and potential of protected areas. Nature, 515, 67–73.
Βασική προϋπόθεση φυσικά είναι να τις έχετε εγκαταστήσει ήδη↩
Πληκτρολογώντας
varsστην κονσόλα σας↩Συγκρίνετε φερειπείν τα δεδομένα που έχουμε για το υψόμετρο με αυτά για οποιαδήποτε άλλη μεταβλητή↩
Δεν λογαριθμούμε τη μεταβλητή
PD. Γιατί;↩Στην προκειμένη περίπτωση αυτό έχει ήδη για εσάς - είναι μια αρκετά επίπονη διαδικασία↩